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Abstract 

The  high  precision  mean  element  (semianalytic)  satellite  theory  developed 
at  Draper  Laboratory  is  more  efficient  than  conventional  numerical  methods  and 
more  accurate  than  the  current  generation  of  analytic  theories.  This  efficiency,  along 
with  its  portability  to  a  variety  of  computing  environments  makes  the  semianalytic 
theory  a  natural  choice  for  maneuver  planning  applications.  These  applications  will 
become  more  important  in  the  future  as  the  capability  of  individual  platforms  to 
maneuver,  the  number  of  platforms  in  space,  and  the  requirements  for  rapid  response 
to  requests  for  data  all  increase.  The  application  of  semianalytic  satellite  theory  to 
an  Earth  Observation  satellite  in  an  orbit  similar  to  that  expected  for  LANDSAT 
6  is  investigated.  Orbit  constraints  such  as  sun  synchronous,  repeat  ground  track, 
frozen  orbit,  and  non-impulsive  maneuver  capabilities  are  included  in  this  analysis. 
Applications  of  maneuver  planning  to  past  and  future  satellite  missions  that  include 
at  least  two  of  the  listed  orbit  constraints  are  discussed. 

Since  atmospheric  drag  is  the  primary  uncertain  disturbing  acceleration  to  the 
nominal  satellite  orbit,  upper  and  lower  limits  of  a  density  confidence  interval  were 
determined.  Two  methods  were  analyzed;  it  was  found  that  using  forecast  and  actual 
solar  flux  and  geomagnetic  activity  data  from  the  years  1986-1990  resulted  in  a  con¬ 
servative  but  realistic  confidence  interval.  The  upper  limit  is  utilized  to  compute  the 
time  of  the  orbit  adjust  burn  and  the  lower  limit  of  the  density  is  used  r  j  calculate 
the  magnitude  of  the  orbit  adjust  burn.  These  limits  are  necessary  so  that  the  ground 
track  boundaries  are  not  exceeded. 
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Chapter  1 
Introduction 

1 . 1  Background 

This  thesis  focuses  on  the  orbital  maneuver  planning  process  for  Earth  Observation 
satellite  missions.  Since  the  launch  of  the  first  Earth  Observation  satellite,  TIROS  -  1, 
in  April  of  1960,  over  135  satellites  have  been  launched  with  specific  missions  to  gain 
data  on  the  Earth  [69] .  Typical  functions  of  Earth  Observation  satellites  encompass 
a  broad  range  including: 


•  Weather 

•  Climate 

•  Meteorology 

•  Land  Survey 

•  Agriculture 

•  Forestry 

•  Disasters  (Alarm,  Relief  Planning) 

•  Hydrology 

•  Oceanology 


A  majority  of  the  measurements  needed  by  these  different  missions  are  common. 
To  attain  these  measurements,  there  are  two  basic  classes  of  instruments,  passive 
and  active.  The  basic  passive  instrument  is  the  radiometer  which  measures  the  re¬ 
flected  radiation  from  the  Earth  within  some  frequency  band  and  polarization.  The 
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microwave  radiometer  detects  sea  ice  conditions  and  cover  in  day  or  night.  It  is 
only  slightly  affected  by  cloud  cover.  The  medium  resolution  visible/infrared  ra¬ 
diometer  can  detect  sea  surface  temperatures,  sea  ice  conditions,  surface  reflectivity, 
ocean  color,  and  classify  vegetation  through  cloud  cover.  The  high  resolution  im¬ 
ager/radiometer,  like  the  thematic  mapper  of  LANDSAT  4/5,  can  determine  changes 
in  vegetation,  ice  motion,  terrain  mapping,  surface  melting  of  snow  and  ice,  and  local 
winds.  In  the  active  instrument  category,  the  most  common  is  the  radar  altimeter, 
which  measures  the  range  to  the  Earth’s  surface.  The  scatterometer  determines  sea 
surface  wind  vectors  from  the  wave  pattern  generated  by  the  wind.  The  synthetic 
aperture  radar,  like  on  SEASAT  and  ERS-1,  is  not  affected  by  lack  of  light  or  cloud 
cover.  It  can  also  penetrate  shallow  ice  or  enow  to  determine  surface  features  below. 
This  precise  instrument  not  only  detects  ice  and  fine  scale  sea  surface  roughness,  but 
also  maps  the  Earth’s  surface  in  all  weather  and  classifies  rock,  soil,  surface  wetness, 
and  vegetation.  One  less  common  active  instrument  is  the  Lidar,  a  laser  altime¬ 
ter  [84],  Some  of  the  instruments  used  frequently  and  their  functions  are  listed  in 
Table  1.1  [3,22,89]. 


Table  1.1.  Earth  Observation  Instrument  Function 


Instrument 

Function 

Altimeter 

•  Surface  topography 

•  Wave  direction 

•  Ocean  wave  height 

Microwave  imager 

•  Wind  speed 

•  Precipitation 

•  Arctic  ice  cover, 
thickness,  and  age 

Microwave  radiometer 

•  Surface  temperatures 

Synthetic  aperture  radar 

•  Image  of  surface 

Scatterometer 

Surface  wind 

16 


The  general  orbital  characteristics  of  these  satellites  repeat  due  to  the  similarity 
in  the  instrument  constraints.  There  are  trade-offs  and  relationships  between  the 
desired  satellite  functions,  the  physical  and  geometrical  constraints  in  the  satellite 
structure,  the  specific  limitations  of  the  vehicle  subsystems,  and  the  properties  of 
orbital  motion  and  satellite  dynamics  [52], 

A  general  Earth  Observation  satellite  will  usually  have  a  mission  lifetime  of  three 
to  five  years.  The  coverage  of  the  Earth  is  usually  global,  but  may  be  constrained 
for  coverage  up  to  a  specified  latitude.  Most  land  observation  satellites  have  coverage 
constraints  between  55°  S  to  75°  N  latitude  due  to  the  distribution  of  the  Earth’s  land 
mass.  An  orbit  of  75°  inclination  covers  the  Earth’s  land  mass  and  leaves  only  the  ex¬ 
treme  north  and  south  latitudes  unobserved.  Other  coverage  constraints  may  include 
a  ground  track  that  repeats  itself.  This  allows  instruments  to  return  to  re-examine  a 
certain  point  of  the  Earth.  A  repeat  ground  track  allows  an  easily  predictable  pat¬ 
tern  of  coverage  and  an  opportunity  to  make  direct  comparisons  between  observations 
taken  at  regular  intervals  for  the  locations  covered.  There  is  a  trade-off  between  the 
frequency  of  re-examination  and  the  degree  of  global  coverage.  The  degree  of  Earth 
coverage  is  also  constrained  by  the  instrument  observation  swath,  or  field  of  view, 
over  the  Earth’s  surface  as  the  satellite  moves  in  its  orbit.  The  swath  coverage  of 
the  Earth  is  determined  by  the  displacement  of  the  swath,  ie.  in  the  ground  track, 
over  time.  The  primary  variable  in  determining  the  swath  width  of  an  instrument  is 
the  altitude  of  the  satellite.  By  raising  or  lowering  the  mean  altitude  of  the  satellite 
orbit,  the  swath  width  respectively  decreases  or  increases.  The  size  of  the  orbit  may 
also  be  determined  by  the  resolution  needed.  If  the  instrument,  itself,  can  not  be 
improved  to  increase  resolution,  the  mean  orbit  altitude  may  be  decreased  to  meet 
the  rosolution  constraints.  But  because  of  the  trade-off  between  mission  lifetime  and 
instrument  resolution,  the  lowering  of  the  mean  orbit  altitude  increases  atmospheric 
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drag  experienced,  thus  decreasing  mission  lifetime  or  increasing  orbital  maintenance 
requirements.  In  addition,  lowering  the  mean  altitude  of  the  satellite  decreases  the 
ground  station  contact  time,  therefore  decreasing  the  amount  of  data  transferred  to  or 
from  the  satellite  [54].  The  conventional  Earth  Observation  satellite  orbit  is  usually 
near  circular  to  ameliorate  altitude  changes  that  may  be  sensed  by  on  board  instru¬ 
ments  over  the  observation  swath.  In  addition,  the  various  aspect  angles  needed  by 
instrument  observations  may  be  a  constraint  on  the  coverage  characteristics  of  an 
Earth  Observation  satellite.  The  instruments  that  attain  these  various  aspect  angles 
must  be  placed  on  board  the  satellite  so  that  the  observations  are  not  contaminated 
by  the  spacecraft  structure.  The  illumination  of  the  Earth  and/or  the  satellite  may  be 
a  deciding  factor  in  determining  the  orbit.  The  existence  of  visible  light  sensors  on  the 
spacecraft  may  cause  a  sun  synchronous  orbit  to  be  desired  for  constant  illumination 
conditions.  As  the  seasons  change,  the  solar  incidence  angle  has  long  periodic  changes 
revealing  landscape  and  terrain  features  on  the  surface  for  a  sun  synchronous  orbit. 
Also,  if  a  satellite’s  power  is  solar  generated,  a  sun  synchronous  orbit  can  minimize 
eclipsing  of  the  solar  wings.  The  number  and  placement  of  sensors  and  experiments 
may  constrain  the  type  of  booster  used,  or  the  booster  weight  and  configuration  limits 
may  constrain  the  instruments  themselves.  The  orbital  maintenance  predicted  for  a 
satellite  in  orbit  may  constrain  the  number  of  solar  cells  on  the  satellite.  Since  larger 
solar  wings  create  a  greater  amount  of  drag,  this  may  decrease  the  lifetime  or  increase 
the  orbital  maintenance  needed.  Another  limit  in  maintenance  burns  is  the  maximum 
velocity  of  the  thrusters  given  by  the  Rocket  Equation: 


where  y  is  the  gravitational  acceleration  of  the  Earth,  Isp  is  the  specific  impulse  of 
the  thruster  fuel,  Winit  is  the  initial  weight  of  the  spacecraft  at  the  beginning  of  the 
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maneuver,  and  Wfinai  is  the  final  weight  of  the  spacecraft  after  the  maneuver  [52], 

To  fully  understand  the  orbital  requirements  of  Earth  Observation  systems,  work¬ 
ing  definitions  of  the  terms  sun  synchronous,  repeat  ground  track,  and  frozen  orbit 
concept  are  needed.  The  detailed  discussion  and  algorithms  associated  with  these 
concepts  will  be  addressed  in  Chapter  4.  Here,  a  brief  preview  of  these  concepts  is 
provided. 

A  sun  synchronous  orbit  keeps  the  orbit  plane  at  a  constant  angle  to  the  sun  by 
using  the  geopotential  perturbations  of  the  Earth  to  move  the  line  of  nodes  westward 
along  the  equator  approximately  one  degree  per  day.  The  line  of  nodes  is  defined  as 
the  intersection  of  the  orbit  plane  and  the  Earth’s  equator.  Through  the  application 
of  the  Earth’s  geopotential  to  Lagrange's  Planetary  Equations,  the  mean  nodal  rate 
is  calculated  by  a  formula  from  Kozai  where  the  elements  are  assumed  to  be  mean 
elements  in  a  secular  sense  [18]: 


dn  3  /  Re  \2  nr 

dt  2  2  ya(l  -  e2) )  Va3C°SZ 


(1.2) 


where  is  the  mean  longitude  of  the  ascending  node,  the  angle  of  the  line  of  nodes 
from  the  vernal  equinox.  The  variable  J2  is  the  harmonic  coefficient  of  the  Earth’s 
oblatcness  equal  to  1.08263xl0-3.  The  mean  equatorial  radius  of  the  Earth,  denoted 
Rn,  is  equal  to  6378.135  km  and  the  gravitational  parameter  of  the  Earth,  p,  is 
398601.2  km3/sec2  [4].  The  exact  values  of  these  parameters  depend  on  the  gravity 
model  employed;  therefore  the  values  stated  here  are  approximate.  The  semi-major 
axis,  a,  the  inclination,  i ,  the  longitude  of  the  ascending  node,  Q,  and  the  eccentricity, 
e,  of  the  orbit  are  variables  in  the  common  Keplerian  orbital  element  set.  For  a  more 
detailed  description  of  orbital  elements,  see  Appendix  A.  By  setting  the  nodal  rate 
equal  to  360°  per  year,  the  solution  for  sun  synchronous  orbital  elements  can  be 
determined  from  Equation  1.2.  To  guarantee  that  the  line  of  nodes  moves  at  the 
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correct  rate,  a  sun  synchronous  Earth  Observation  satellite  is  usually  in  a  slightly 
retrograde  orbit  with  an  altitude  of  approximately  200  to  1000  nautical  miles  (370  to 
1852  km). 

A  satellite  ground  track  is  defined  as  the  locus  of  points  traced  out  on  the  Earth’s 
surface  directly  beneath  the  spacecraft  orbit.  This  ground  track  does  not  normally 
repeat  itself,  ie.  retrace  a  previous  orbit’s  ground  track,  by  reason  of  the  Earth 
rotating  underneath  the  satellite  orbit  while  the  line  of  nodes  rotates  around  the 
Earth’s  equator  due  to  perturbations.  A  commensurability  condition  for  a  repeat 
ground  track  may  be  determined  by  using  an  equation  from  Baxter  [5]: 

2irN  =  M(ue  -  f l)PN  (1.3) 

where  N  is  the  number  of  nodal  days  until  the  ground  track  repeats,  M  is  the  number 
of  orbits  in  the  repeat  cycle,  u ;e  is  the  approximate  constant  rotation  rate  of  the 
Earth,  fi  is  the  rotation  rate  of  the  line  of  nodes,  and  Pjv  is  the  nodal  period  of  the 
satellite.  A  nodal  day  is  defined  as  the  time  for  the  Earth  to  rotate  360°  with  respect 
to  the  line  of  nodes.  A  solar  day  is  defined  as  the  time  for  the  Earth  to  rotate  360° 
with  respect  to  the  Earth-Sun  line.  If  the  satellite  orbit  is  sun  synchronous,  a  solar 
day  is  approximately  the  same  as  a  nodal  day  since  the  line  of  nodes  is  at  a  constant 
angle  to  the  sun.  The  nodal  period  depends  mainly  on  the  Keplerian  period,  orbit 
inclination,  and  to  a  lesser  extent  the  orbit  eccentricity.  If  the  satellite  is  near  circular, 
ie.  eccentricity  is  approximately  zero,  the  equation  for  the  nodal  period  is  [78,18]: 

PN  =  Pk  1  -  ~.h  )  (6  -  5 sin2  i)  +  0(e2)  (1.4) 

where  f\  is  defined  as  the  Keplerian  (two  body)  orbital  period  [4]: 
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where  a  is  the  semi-major  axis  and  i  is  the  inclination  and  both  values  are  taken  at  the 
node  crossing.  As  previously  stated  in  Equation  1.2,  the  variable  J2  is  the  harmonic 
coefficient  of  the  Earth’s  oblateness  equal  to  1.08263xl0~3,  the  mean  equatorial  ra¬ 
dius  of  the  Earth,  denoted  Re,  is  equal  to  6378.135  km,  and  the  Earth’s  gravitational 
parameter,  p,  is  398601.2  km3/sec2.  Again,  these  values  are  approximate  and  depen¬ 
dent  on  the  actual  gravity  model  used.  A  repeat  ground  track  orbit  is  advantageous 
because  it  allows  observations  of  areas  of  interest  to  be  regularly  scheduled. 

The  frozen  orbit  concept  sets  the  argument  of  perigee  to  ±90°  and  the  eccentric¬ 
ity  to  a  determined  small  nominal  value.  This  configuration  causes  the  J2  and  odd 
numbered  harmonics  of  the  Earth’s  gravitational  field  to  cancel  one  another  causing 
the  argument  of  perigee  to  oscillate  around  the  nominal  value  if  it  is  perturbed.  By 
setting  the  argument  of  perigee  to  ±90°,  the  eccentricity  oscillates  around  its  set  nom¬ 
inal  value  since  the  lower  order  geopotential  perturbations  in  the  eccentricity  depend 
on  the  cosine  of  the  argument  of  perigee.  The  frozen  orbit  concept  provides  passive 
eccentricity  and  argument  of  perigee  control.  In  addition,  it  minimizes  spacecraft 
altitude  variations  at  all  latitudes  since  the  perigee  location  is  nearly  constant  [61]. 

The  optimal  compromise  between  these  parameters  to  accomplish  the  Earth  Ob¬ 
servation  mission  is  created  using  two  body  orbital  mechanics  plus  the  oblateness  of 
the  Earth  to  identify  nominal  mission  orbits.  Then,  other  orbital  perturbations  are 
added  in  to  determine  if  the  candidate  orbit  will  realistically  succeed.  By  determin¬ 
ing  how  this  orbit  will  change  over  time,  the  maintenance  maneuvers  needed  may  be 
accurately  planned  before  launch  to  determine  the  fuel  and  maneuver  schedule  for 
a  successful  mission.  This  candidate  orbit  must  be  feasible  in  terms  of  the  mainte¬ 
nance  required  and  its  accessibility  for  the  launch  vehicle  expected.  The  maintenance 
needed  must  not  be  too  expensive  in  terms  of  fuel  weight  for  the  desired  mission  in 
addition  to  not  interfering  with  the  observations  by  the  instruments  on  board  the 
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spacecraft.  If  any  one  of  these  constraints  is  violated,  a  new  nominal  orbit  must  be 
determined  and  the  orbit  analysis  cycle  repeated. 


1.2  Maneuver  Planning  for  Single  Satellites 

There  are  several  past  and  future  Earth  Observation  satellite  applications  that 
jointly  utilize  at  least  two  of  the  following  concepts:  sun  synchronous,  repeat  ground 
track,  and  frozen  orbit.  Past  missions  include  the  1978  precursor  SEASAT,  the 
LANDSAT  4/5  missions  in  1982  and  1984  respectively,  the  NROSS  program  of  1986 
that  was  killed  prematurely,  and  GEOSAT  from  1986  to  1990.  Even  though  NROSS 
was  not  launched,  it  motivated  several  important  analytical  efforts  and  is  considered 
here.  In  addition,  there  are  several  future  programs  that  will  follow  similar  orbits. 
Among  these  are  ERS-1,  LANDSAT  6,  TOPEX,  and  EOS.  A  more  specific  outline 
of  these  missions  follows. 

1.2.1  SEASAT 

After  SEASAT’s  launch  on  June  26,  1978,  NASA  proved  that  oceanographic 
parameters  could  be  measured  from  a  satellite  observatory.  Its  accurate  observations 
with  a  radar  altimeter  and  microwave  sensors  included  ocean  surface  temperatures, 
winds,  waves,  sea  ice,  currents,  and  atmospheric  water  content  [21,71,89].  SEASAT 
was  designed  for  a  one  year  life  time  with  the  first  two  months  used  to  insert  the 
satellite  in  orbit  and  calibrate  the  instruments.  This  initial  orbit  was  an  exact  three 
day  repeat  ground  track,  renewing  every  43  orbit  revolutions.  This  ground  track 
passed  over  the  laser  ranging  site  on  the  island  of  Bermuda  to  calibrate  the  altimeter. 
The  rapid  global  sampling  achieved  in  this  orbit  also  was  advantageous  for  calibrating 
the  other  instruments  on  board  SEASAT.  In  mid-August  to  accomplish  the  mission 
objectives,  the  orbit  was  altered  to  a  near  ‘Cambridge’  orbit  with  a  repeat  cycle  that 
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lasted  17  nodal  days  or  244  orbit  revolutions  [69,89].  Due  to  a  malfunction  in  the 
power  subsystem,  SEAS  AT  failed  on  October  9,  1978  [3].  The  mean  orbital  elements 
of  SEASAT  while  it  was  successfully  in  orbit  are  stated  in  Table  1.2  [21,89]. 


Table  1.2.  Orbital  Elements  for  SEASAT 


Mean  Keplerian  Element 

SEASAT  Injection 

SEASAT  Mission 

semi-major  axis,  a 

7168.3  km 

7173.6  km 

eccentricity,  e 

«  0.0008 

«  0.001045786 

inclination,  i 

108° 

108° 

argument  of  perigee,  u 

90° 

90° 

These  orbital  elements  supported  the  two  separate  repeat  ground  tracks  in  a  frozen 
orbit  condition.  The  mission  orbit  mapped  95%  of  the  Earth’s  global  oceans  every 
thirty-six  hours  [69].  The  results  achieved  before  SEASAT  failed  demonstrated  the 
feasibility  and  necessity  of  using  satellite  observations  to  obtain  global  information. 

1.2.2  LANDSAT  4/5 

The  Earth  observation  satellites,  LANDSAT  4  and  LANDSAT  5,  were  both 
launched  from  the  Western  Test  range  on  Delta  3920  launch  vehicles  [42].  Both 
LANDSAT  4  and  5  had  two  instruments  to  make  observations  of  the  Earth.  The 
multi-spectral  scanning  radiometer  and  the  high  resolution  thematic  mapper  were 
used  to  take  surface  pictures  of  the  Earth  [69],  There  were  many  applications  for 
the  observations  from  the  LANDSAT  program:  agribusiness,  geology,  forestry,  disas¬ 
ter  assessment/engineer  planning,  hydrology,  land  use  and  regional  planning,  range 
management,  and  cartography  [20],  After  LANDSAT  4  encountered  difficulties  in  re¬ 
turning  data  from  its  thematic  mapper,  its  multi-spectral  scanner  was  used  to  obtain 
data  for  foreign  ground  stations  while  LANDSAT  5  became  the  global  provider  for 
thematic  mapping  data  [69]. 
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LANDSAT  4  was  put  into  orbit  first  in  July  of  1982,  and  then  LANDSAT  5  fol¬ 
lowed  on  March  1,  1984.  LANDSAT  5  was  put  in  the  same  orbit  and  phased  180°  with 
LANDSAT  4  so  that  their  period  of  world  coverage  was  halved.  Together  LANDSAT 
4  and  5  form  a  two  satellite  constellation  but  each  is  treatr  ’  as  an  individual  satel¬ 
lite.  This  is  done  by  maintaining  each  satellite  in  a  sun  synchronous,  repeat  ground 
track,  frozen  orbit  [42,69].  The  orbit  inclination  maintenance  preserves  the  sun  syn¬ 
chronous  orbit  which  maintains  the  nodal  period  needed  for  a  repeat  ground  track. 
The  repeat  ground  track  maintenance,  within  ±10  km,  preserves  the  common  orbit 
that  the  satellites  share  and  the  phasing  between  the  two  satellites.  The  frozen  orbit 
concept  allows  for  passive  maintenance  of  the  eccentricity  and  argument  of  perigee. 
The  frozen  orbit  concept  also  favors  ideal  conditions  for  the  instruments  on  board 
since  the  altitude  variations  are  minimized  [49].  Both  satellites  are  in  a  16  day  re¬ 
peat  ground  track  cycle.  Each  satellite  achieves  giobal  coverage  in  one  16  day  repeat 
cycle,  or  every  233  revolutions  of  the  satellite  in  orbit.  Since  LANDSAT  5  is  phased 
8  days  from  LANDSAT  4,  there  is  global  coverage  every  8  days  if  both  satellites  are 
operational.  The  repeat  ground  track,  sun  synchronous,  and  frozen  orbit  is  achieved 
by  the  orbital  elements  shown  in  Table  1.3  [42]. 


Table  1.3.  Orbital  Elements  for  LANDSAT  4/5 


Kcplerian  Element 

LANDSAT  4/5  Value 

semi-ma  jor  axis,  a 

7077.8  km 

mean  eccentricity,  e 

^  0.0012 

inclination,  i 

98.2° 

argument  of  perigee, 

90.0° 

equatorial  crossing  time 

0930-1000  hours  local  time 

Orbital  perturbations  cause  these  nominal  orbital  elements  to  vary  from  their  ini¬ 
tial  values.  'Io  predict  the  frequency  of  the  altitude  and  inclination  adjust  burns, 
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LANDS  AT  mission  analysis  used  a  Cowell  propagator  to  generate  an  accurate  pre¬ 
diction  of  the  LANDS  AT  4/5  orbit  [42],  Atmospheric  drag  is  the  primary  perturbing 
force  for  a  satellite  in  a  low  altitude  orbit.  Since  the  in-plane  force  of  drag  causes 
the  satellite  to  lose  energy,  the  satellite  loses  orbital  altitude  and  the  satellite  period 
decreases,  see  Equation  1.5.  This  can  be  viewed  as  the  satellite  arriving  at  its  equator 
crossing  too  early.  This  does  not  allow  the  Earth  enough  time  to  rotate  sufficiently  for 
the  ground  track  to  exactly  repeat,  thus  the  ground  track  appears  to  drift  eastward. 
To  keep  the  ground  track  within  its  required  boundaries  of  ±  10  km,  an  altitude  adjust 
burn  must  be  done  approximately  every  50  days  for  a  median  value  of  atmospheric 
drag.  This  median  atmosphere  was  based  on  the  average  of  the  estimated  extremes  of 
the  solar  flux  values  during  the  mission.  A  graphical  view  of  this  maneuver  estimation 
can  be  viewed  in  the  LANDSAT-D  Orbit  Adjust  Criteria  [49].  The  atmospheric  drag 
is  determined  mathematically  through,  an  atmospheric  model  and  the  prediction  of 
the  solar  flux  values  in  the  future.  These  solar  flux  values  are  difficult  to  forecast 
and  are  discussed  fully  in  Chapter  5.  If  the  solar  flux  is  less  than  predicted  after 
the  maneuver,  the  ground  track  may  drift  beyond  the  western  boundary  requiring  a 
retrograde  maneuver  (a  180°  spacecraft  yaw).  Therefore,  the  altitude  adjust  burns 
must  be  conservative  in  estimating  the  burn  magnitude  [42]. 


The  inclination  drift  is  primarily  due  to  solar  gravitational  forces  acting  on  the 
satellite.  These  forces  are  not  in  the  orbital  plane  and  cause  the  orbit  to  become  more 
polar.  This  change  in  inclination  affects  the  sun  synchronous  condition  of  the  orbit, 
thus  necessitating  adjustment  burns  [49].  The  orbital  inclination  after  orbit  injection 
was  set  to  a  value  such  that  it  was  18  months  after  launch  that  the  first  inclination 
adjustment  burn  was  needed.  After  that,  an  inclination  burn  was  anticipated  to  be 
needed  every  8  months  [42] . 
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1.2.3  Navy  Remote  Ocean  Sensing  System  (NROSS) 


The  NROSS  mission,  a  derivative  of  the  SEASAT  mission,  was  to  map  sea 
surface  topography  by  measuring  mesoscale  variations  in  the  sea  surface  height.  Sea 
surface  topography  is  computed  as  the  residual  of  the  reference  sea  geoid  and  the 
ocean  surface  height  detected  by  the  on  board  altimeter.  This  altimeter  is  identical 
to  the  altimeter  presently  on  board  the  GEOSAT  spacecraft  [3],  The  sea  geoid  is 
the  expected  ocean  surface  affected  only  by  the  Earth’s  gravitational  potential.  This 
includes  the  reference  ellipsoidal  earth  plus  the  sea  height  variation  due  to  the  tides. 
An  illustration  of  the  definition  of  sea  surface  topography  can  be  seen  in  Figure  1.1. 
The  value,  H,  is  the  observed  altimeter  height  of  NROSS  while  alt  is  the  radial 
altitude  of  the  satellite  independently  determined  from  the  altimeter  height.  The 
difference  between  these  two  values  gives  the  ocean  topography  [68].  Since  the  ocean 


ellipsoid 

Figure  1.1.  Illustration  of  Sea  Surface  Topography 
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geoid  is  not  accurately  known  over  most  regions,  NROSS  was  to  primarily  track 
changes  in  rather  than  the  absolute  value  of  the  ocean  topography  [68]. 

The  NROSS  spacecraft  was  a  scaled  down  version  of  the  cancelled  National  Oceanic 
Satellite  System  (NOSS)  from  1981.  The  NOSS  participants,  ie.  NASA,  Navy,  and 
National  Oceanic  and  Atmospheric  Administration  (NOAA),  were  joined  again  un¬ 
der  the  leadership  of  the  Navy  to  form  the  NROSS  program  in  1985.  The  Air  Force 
also  joined  in  to  provide  ground  station  and  tracking  facilities  [71].  In  one  view, 
the  NROSS  spacecraft  was  to  be  a  NOAA-D  platform  bus  with  Defense  Meteorlogi- 
cal  Satellite  Program  (DMSP)  communications  and  telemetry  processing  subsystems. 
The  NROSS  instruments  were  to  be  added  on  this  bus  to  accomplish  the  mission  ob¬ 
jectives  [22].  The  instruments  scheduled  to  be  on  board  NROSS  were  the  altimeter 
previously  mentioned,  a  modified  scatterometer  from  SEASAT,  a  new  low-frequency 
radiometer,  and  a  special  sensor  microwave/imager  (SSM/I)  from  DMSP  rather  than 
the  synthetic  aperture  radar  from  SEASAT  [3,69]. 

The  optimal  orbit  called  for  a  sun  synchronous  and  stringent  repeat  ground  track 
orbit  [18].  The  repeat  ground  track  lessened  altimeter  bias  from  ocean  geoid  height 
at  points  where  the  current  ground  track  crossed  over  a  previous  ground  track  [68]. 
Originally,  a  0.001  maximum  mean  eccentricity  constraint  was  specified.  This  near 
circular  orbit  was  to  be  maintained  throughout  the  mission  to  avoid  unnecessary  al¬ 
titude  variations  in  instrument  observations  of  the  sea  surface  height.  The  on  board 
altimeter  experiment  needed  the  support  of  a  rigorous  repeat  ground  track  and  the 
scatterometer  required  minimal  altitude  variations  to  maintain  stable  locations  for 
the  sea  surface  cells.  Using  the  frozen  orbit  concept  increased  the  nominal  mean 
eccentricity  to  0.00115  but  minimized  altitude  variations  throughout  all  latitudes  of 
the  orbit.  The  frozen  orbit  concept  also  had  an  additional  advantage  by  decreas¬ 
ing  the  amount  of  orbit  maintenance  through  passive  control  of  the  eccentricity  and 
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argument  of  perigee.  Hence,  a  frozen  orbit  geometry  was  adopted,  since  the  frozen 
orbit  advantages  pertaining  to  decreased  orbit  maintenance  and  increased  instrument 
performance  outweighed  the  small  increase  in  maximum  mean  eccentricity.  The  or¬ 
bital  elements  prescribed  for  the  NROSS  mission  are  summarized  in  Table  1.4.  These 


Table  1.4.  Orbital  Elements  for  NROSS 


Keplerian  Element 

NROSS  Value 

semi-major  axis,  a 

7198.65  km 

mean  eccentricity,  e 

0.00115 

inclination,  i 

98.697° 

argument  of  perigee,  u 

90.0° 

equatorial  crossing  time 

0530  hours  local  time 

elements  insured  a  frozen,  sun  synchronous  orbit  and  a  repeat  cycle  of  19  nodal  days, 
or  270  orbit  revolutions  [61]. 

NROSS  was  to  be  launched  on  a  Titan  II  launch  vehicle  near  the  maximum  of  the 
11  year  solar  cycle,  which  was  expected  in  1990.  These  maximum  solar  flux  values 
were  expected  to  create  high  levels  of  atmospheric  drag  during  the  three  to  four  year 
mission  lifetime  of  NROSS  [34,3,71].  In-plane,  along  track  burns  maintain  both  the 
semi-major  axis  and  the  orbital  eccentricity  and  are  designed  to  simultaneously  pre¬ 
serve  both  the  frozen  orbit  geometry  and  meet  the  repeat  ground  track  requirement. 
Because  of  high  estimated  solar  activity,  the  semi-major  axis  decay  due  to  atmospheric 
drag  was  the  expected  dominant  cause  of  ground  track  drift.  Because  of  this  drift, 
altitude  adjust  burns  were  predicted  to  take  place  every  10-25  days  [33].  In  addition, 
inclination  drift  causes  ground  track  error  at  extreme  latitudes.  It  was  determined  in 
the  NROSS  orbit  adjust  strategy  that  a  bound  on  the  inclination  of  ±0.009°  would 
maintain  the  ground  track  at  the  high  latitudes  [34].  Additional  ground  track  drift 
was  attributed  to  changes  in  the  nodal  rate  and  the  mean  eccentricity. 
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The  NROSS  project  had  its  procurement  cancelled  in  December  of  1986,  and  the 
entire  project  cancelled  in  1987  [61].  The  mission  objectives  of  NROSS  have  been  split 
up  and  partitioned  to  other  satellite  systems,  ie.  scatterometer  on  EOS,  altimeter  on 
SPIN  satellite,  and  an  altimeter  and  radar  scatterometer  on  DMSP  Block  6  [73]. 

1.2.4  GEOSAT 

The  United  States  Navy  GEOSAT  spacecraft  was  launched  March  12,  1985, 
on  an  Atlas  E  launch  vehicle  to  complete  geoid  mapping  initialized  by  SEASAT  in 
1978  [53].  After  this  mission  was  accomplished,  GEOSAT  performed  a  transition  burn 
between  October  1,  1986  and  November  6,  1986,  from  a  geodetic  to  an  oceanographic 
phase  of  the  mission.  This  second  phase  is  referred  to  as  the  GEOSAT  Exact  Repeat 
Mission  (ERM)  [30].  The  only  instrument  on  board  was  a  highly  redundant  altimeter. 
Although  this  instrument  had  a  narrow  swath  width  causing  the  data  collected  to  be 
exiguous,  the  altimeter  observed  ocean  wave  height  and  direction  and  tracked  ocean 
surface  topography.  The  wind  speed  and  direction  could  be  determined  from  ground 
processing  of  the  received  back  scatter  cross  section  [3,53].  The  altimeter  data  was 
also  used  to  create  a  dense  grid  of  observations  required  to  improve  the  model  of  the 
Earth’s  gravitational  potential  [69].  The  GEOSAT  ERM  ground  track  repeated  every 
244  orbits  or  every  17  nodal  days  and  was  maintained  until  1990  within  ±  1  km.  The 
research  initialized  by  GEOSAT  will  be  followed  on  by  SALT  in  1991  [30]. 

GEOSAT  measured  the  time  dependent  oceanographic  features,  ie.  rings  and 
eddies,  and  the  time  independent  features  of  circulation,  ie.  gulf  stream,  in  addition 
to  the  sea  geoid  height.  It  was  necessary  to  implement  a  stringent  repeat  cycle 
specification  to  minimize  the  effects  of  sea  geoid  uncertainty  on  the  determination  of 
ocean  variability,  by  allowing  the  sea  geoid  to  be  averaged  out  of  the  satellite  altimeter 
data  [53].  GEOSAT  also  utilized  the  frozen  orbit  concept  resulting  in  a  constant 
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altitude  history  from  orbit  to  orbit  since  the  argument  of  perigee  and  eccentricity 
were  stable  [79], 

To  maintain  the  repeatable,  frozen  orbit  condition,  the  ERM  orbital  elements  are 
summarized  in  Table  1.5.  These  elements  gave  GEOSAT  a  nodal  rate  of  2. 0517209° /day 

Table  1.5.  Orbital  Elements  for  GEOSAT  ERM 


Keplerian  Element 

GEOSAT  ERM  Value 

semi-major  axis,  a 

7173.6  km 

mean  eccentricity,  e 

0.0007971 

inclination,  i 

108.04397° 

argument  of  perigee,  cj 

91.493632° 

longitude  of  the  ascending  node,  Q 

49.788597° 

and  a  nodal  period  of  6037.5601  seconds  [79]. 

GEOSAT  maneuver  research  determined  that  a  change  in  the  semi-major  axis 
is  most  efficient  with  an  along  track  burn  and  can  be  changed  independent  of  the 
other  Keplerian  orbital  elements  if  equal  burns  are  done  at  any  true  anomaly,  /,  and 
/  ±  180°.  A  change  in  the  eccentricity  of  the  orbit  is  also  most  efficient  with  an 
along  track  burn  and  can  be  changed  independent  of  the  other  elements  with  equal 
and  opposite  burns  at  /  =  0°  and  /  =  180°.  A  change  in  the  argument  of  perigee 
is  also  most  efficient  with  a  burn  in  the  along  track  direction  and  independent  of 
the  other  elements  with  equal  and  opposite  burns  at  true  anomalies  of  /  =  ±90°. 
A  simultaneous  change  in  the  argument  of  perigee  and  eccentricity  can  be  achieved 
independent  of  a  change  in  the  semi-major  axis  with  along  track  burns  of  equal  and 
opposite  magnitudes  at  true  anomalies,  /  and  /  +  180°,  where  /  is  determined  by  the 
equation: 

/  =  -u  +  arctan  [(esinw  —  e0)/ecosu;]  (1.6) 

where  e  is  the  actual  eccentricity  of  the  orbit,  c0  is  the  nominal  frozen  orbit  eccentric- 


ity,  and  u  is  the  argument  of  perigee  [79].  These  general  independent  element  adjust 
burns  are  summarized  in  Table  1.6. 

Table  1.6.  Independent  Element  Adjust  Burns 


True  Anomaly,  / 

Burn  Direction 

a 

fi\ h  =  fi  ±  180° 

Both  posigrade 

e 

O 

O 

oo 

r— H 

i! 

O  . 

O  | 

ii 

5 

one  posigrade,  one  retrograde 

U) 

/i  =  90°;  /2  =  -90° 

one  posigrade,  one  retrograde 

Simulations  and  real  time  data  show  that  the  ground  track  drift  is  sensitive  to 
errors  in  nodal  period,  node  rate,  and  atmospheric  drag,  with  drag  being  the  least 
predictable  of  the  three.  The  average  decrease  in  the  semi-major  axis  during  the 
initial  geodetic  mission  was  approximately  0.5  m/day,  caused  by  an  Aio.r  solar  flux 
between  70  and  75.  Assuming  that  the  solar  flux  would  be  approximately  the  same 
for  the  ERM,  the  time  of  arrival  at  the  ascending  node,  after  N  revolutions,  is  reduced 
in  seconds  per  day  by: 

A tN  =  2.2  x  10"5iV2  (1.7) 

This  equation  is  derived  from  taking  the  partial  derivative  of  the  nodal  period  with  re¬ 
spect  to  the  semi-major  axis  and  assuming  the  ERM  semi-major  axis  to  be  7167.4  km. 
To  maintain  the  ERM  within  one  kilometer,  the  nodal  crossing  needs  to  be  maintained 
within  ±2.15  sec  of  the  ERM  nodal  crossing  time.  This  means  a  theoretical  22  days 
until  the  first  burn  is  needed  at  the  eastern  boundary.  The  orbit  is  then  adjusted 
by  overshooting  the  semi-major  axis  needed  for  the  exact  repeat  period  and  allowing 
the  ground  track  to  drift  west.  Orbital  perturbations,  mainly  drag,  will  decrease  the 
semi-major  axis  and  force  the  ground  track  through  the  exact  repeat  period.  An¬ 
other  maneuver  will  be  performed  before  the  ground  track  violates  the  eastern  1  km 
boundary.  Mathematically,  this  would  lead  to  approximately  50  days  between  follow- 
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ing  orbit  adjust  maneuvers,  even  with  a  10%  error  in  drag  prediction  [9j.  The  actual 
maintenance  maneuvers  for  the  semi-major  axis  were  done  approximately  every  30 
days,  since  the  solar  flux  varied  more  than  expected  [79]. 

1.2.5  ERS-1 

The  European  Space  Agency  (ESA)  will  launch  its  first  remote  sensing  spacecraft, 
ERS-1,  on  May  3,  1991.  Its  mission  is  to  observe  coastal  oceans  including  ice  forma¬ 
tions,  while  measuring  ocean  and  wind  data  to  improve  global  weather  information, 
thus  making  it  an  operational  successor  to  SEASAT  and  NROSS.  The  information 
received  from  this  satellite  will  complement  the  optical  data  from  the  LANDSAT  and 
French  SPOT  systems.  It  is  expected  that  ERS-1  will  have  a  gradual  transfer  from 
an  initial  experimental  mode  to  an  operational  mission  to  attain  usable  Earth  obser¬ 
vations  [69,10].  The  instruments  that  ERS-1  will  carry  to  accomplish  this  mission 
arc  listed  in  Table  1.7.  The  functions  of  these  instruments  can  be  found  in  Table  1.1. 

Table  1.7.  ERS-1  Instrumentation 

•  Wind  scatterometer 

•  Radar  altimeter 

•  Synthetic  aperture  radar 

•  Laser  retro-reflector  array 

•  Along-track  scanning  radiometer 

•  Precise  range  and  range  rate  exp. 


The  ERS-1  synthetic  aperture  radar  (SAR)  has  two  operational  modes.  The  wide 
swath  width  of  the  Image  Mode  observes  all  weather  images  over  oceans,  polar  re¬ 
gions,  and  land.  The  Wave  Mode  of  the  SAR  creates  5  km  x  5  km  images  at  regular 
intervals  to  determine  the  length  and  direction  of  ocean  waves.  The  European  Space 
Operations  Center  (ESOC)  will  determine  orbit  data  within  60  m  along  track  through 
accurate  modelling  of  atmospheric  drag  and  the  Earth’s  geopotential.  The  data  from 
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the  laser  and  precise  range  and  range  rate  experiment  (PRARE)  will  ameliorate  this 
orbit  determination  and  increase  autonavigation  [41,2,10].  The  PRARE,  sponsored 
by  the  German  Federal  Ministry  for  Research  and  Technology,  will  have  its  perfor¬ 
mance  tested  on  the  ERS-1  mission  for  future  utilization  on  other  space  missions. 
The  PRARE’s  space  segment,  which  includes  its  own  data  transmission  and  memory 
for  a  more  generic  system,  sends  a  signal  down  to  an  unmanned  tracking  station.  The 
signal  is  radiated  at  frequencies  of  8  and  2  GHz.  The  ground  station  works  as  a  re¬ 
generative  transponder  for  the  8  GHz  signal.  The  primary  measurement  is  performed 
bv  the  space  segment  at  reception  of  the  returning  7  GHz  signal.  The  accuracy  is 
expected  to  be  within  0.1  -  0.2  m  for  range  and  within  0.3  mm/sec  for  range  rate. 
Since  the  2  GHz  signal  is  affected  more  by  the  ionosphere,  there  is  a  time  difference 
(with  the  8  GHz  signal)  measured  at  the  ground  station  which  is  also  returned  to  the 
satellite.  At  some  central  ground  station,  the  time  difference  is  processed  to  estimate 
the  ionospheric  correction  to  the  range  and  range-rate  data  to  be  used  in  the  orbit 
determination.  The  PRARE  system  will  function  in  all  weather  since  the  signals  are 
in  the  microwave  band  of  frequencies  [40]. 

All  of  the  instruments  listed  in  Table  1.7  need  to  observe  the  fast  changing  features 
over  the  globe  in  constant  geometric  and  local  time  conditions  to  accomplish  their 
mission.  This  leads  to  a  trade-off  between  frequency  of  re-observation  of  covered  areas 
and  fraction  of  global  coverage.  The  orbital  elements  chosen  to  meet  these  needs 
are  listed  in  Table  1.8.  This  orbit  is  sun  synchronous  to  provide  regular  observation 
times  and  constant  lighting,  especially  for  the  along-track  scanning  radiometer.  Three 
separate  repeat  ground  tracks  will  be  supported  by  these  orbital  elements  during  the 
2-3  year  life  time  of  ERS-1.  The  first  repeat  cycle  will  be  every  three  days,  the 
second  repeat  cycle  will  last  35  days,  and  the  third  will  last  176  days.  All  of  the 
ground  tracks  will  be  maintained  within  ±1  km.  The  altitude  adjusts  needed  to 
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Table  1.8.  Orbital  Elements  for  ERS-1 


Mean  Keplerian  Element 

ERS-1  Value 

average  semi-major  axis,  a 

7153.1439  km 

eccentricity,  e 

*  0.001166 

inclination,  i 

98.52146° 

argument  of  perigee,  u 

90.0° 

descending  equatorial  crossing  time 

1030  hrs  local 

maintain  these  ground  tracks  will  be  provided  by  thrusters  in  the  posigrade  direction. 
ERS-1  will  be  yawed  90°  to  provide  an  out  of  plane  thrust  for  periodic  inclination 
adjusts.  The  central  ground  station  at  Salmijaervi,  Sweden,  will  link  up  with  ERS-1 
approximately  10  times  per  day.  At  least  once  a  day,  the  ground  station  will  verify 
satellite  operations  scheduled  for  the  next  24  hrs,  including  necessary  orbit  adjusts. 
The  repeat  ground  track  coverage  attained  depends  upon  both  the  field  of  view  of  the 
instruments  and  the  mean  spacing  between  adjacent,  not  successive,  ground  tracks. 
For  the  ERS-1  three  day  repeat  ground  track  orbit,  the  distance  between  adjacent 
ground  tracks  is  predicted  to  be  932  km.  The  argument  of  perigee  was  set  to  90.0° 
and  the  eccentricity  calculated  to  attain  a  frozen  orbit  to  minimize  altitude  variations 
and  provide  passive  eccentricity  and  argument  of  perigee  control  [41].  It  is  expected 
if  solar  activity  is  low,  F10.7  =  70  and  Ap  —  8,  maneuvers  to  maintain  these  elements 
will  take  place  approximately  every  7  weeks.  If  the  solar  activity  is  high,  F10. 7  =  210 
and  Ap  =  10,  maneuvers  are  expected  every  2-4  weeks  [89,35,2]. 

The  follow  on,  ERS-2,  will  be  launched  in  1994  by  an  Ariane  4  launch  vehicle. 
This  satellite  will  continue  the  data  observed  by  ERS-1  with  the  addition  of  a  global 
ozone  monitoring  experiment  to  examine  the  ozone  problem  [41]. 
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1.2.6  LANDSAT  6 


LANDSAT  6  is  an  advanced  Earth  Observation  satellite  to  be  placed  in  orbit  in 
May,  1992,  with  a  Titan  II  launch  vehicle  [65].  Like  LANDSAT  4/5,  LANDSAT  6 
will  be  in  a  sun  synchronous,  repeat  ground  track,  frozen  orbit.  The  orbital  elements 
that  accomplish  these  orbit  constraints,  listed  in  Table  1.9,  are  very  similar  to  the 
LANDSAT  4/5  elements,  listed  in  Table  1.3  [70,62]. 


Table  1.9.  Orbital  Elements  for  LANDSAT  6 


Keplerian  Element 

LANDSAT  6  Value 

semi- major  axis,  a 

7077.8  km 

mean  eccentricity,  e 

«  0.0013 

inclination,  i 

98.2° 

argument  of  perigee,  uj 

90.0° 

equatorial  crossing  time 

0930-0945  hours  local  time 

As  in  the  previous  LANDSAT  4/5  satellites,  the  LANDSAT  6  ground  track  will 
repeat  every  233  orbits  or  every  sixteen  days.  Thus,  adjacent,  not  necessarily  succes¬ 
sive,  ground  tracks  will  be  separated  by  approximately  1.54506°.  This  ground  track 
repeat  cycle  will  be  maintained  within  ±5  km  during  the  desired  5  year  lifetime. 
Since  the  time  of  launch  places  the  orbit  near  the  peak  of  the  solar  cycle,  atmospheric 
drag  will  be  the  primary  perturbation  to  counteract  with  altitude  adjusts  to  maintain 
the  repeat  ground  track.  The  needed  change  in  eccentricity  to  maintain  the  frozen 
orbit  condition  will  most  often  be  corrected  at  the  same  time  as  the  altitude  adjust. 
If  the  change  in  eccentricity  is  too  large  to  be  completely  reset  during  the  altitude 
adjust,  the  eccentricity  adjust  will  be  of  the  greatest  degree  possible  given  only  posi- 
grade  burns  are  allowed.  As  previously  observed  in  Table  1.6,  an  independent  change 
in  eccentricity  is  accomplished  with  burns  in  equal  and  opposite  directions.  Since 
the  LANDSAT  6  spacecraft  bus  will  only  have  thrusters  in  the  posigradc  direction, 
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the  desired  retrograde  burn  to  independently  change  the  eccentricity  would  require  a 
costly  180°  yaw  or  pitch  to  position  the  thrusters  in  the  correct  burn  direction.  The 
yaw  maneuver  is  preferred  so  that  the  satellite  will  not  flip  over  180°  as  would  occur 
in  the  undesirable,  but  possible,  pitch  maneuver.  In  addition,  an  out  of  plane  burn 
will  require  an  approximately  90°  yaw  to  correctly  position  the  thrusters  to  adjust 
the  inclination.  The  inclination  will  need  to  be  adjusted  periodically  to  maintain  the 
sun  synchronous  condition  and  the  repeat  ground  track  at  high  latitudes.  Since  the 
primary  perturbation  that  changes  the  orbit  inclination  is  due  to  solar  gravity  and  is 
nearly  constant  for  a  frozen,  sun  synchronous  orbit  [65,70,62],  the  inclination  adjusts 
are  simple  to  predict.  The  inclination  adjust  needed  will  be  similar  to  the  LANDSAT 
5  adjusts  at  a  A i  less  than  0.03°. 

Since  the  government  is  now  commercializing  the  LANDSAT  mission  with  gov¬ 
ernment  funding,  LANDSAT  6  is  being  designed  and  built  by  General  Electric  Astro- 
Space  Division  for  the  Earth  Observation  Satellite  (EOSAT)  Company.  The  National 
Oceanographic  and  Atmospheric  Administration  (NOAA)  oversees  the  LANDSAT  6 
project  since  they  managed  the  LANDSAT  4  and  LANDSAT  5  projects  for  the  US 
government  [58,65]. 

1.2.7  Ocean  Topography  Experiment  (TOPEX) 

The  joint  US/French  TOPEX/POSEIDON  will  be  launched  in  June,  1992,  by  the 
National  Aeronautics  and  Space  Administration  on  a  French  Ariane  4  launch  vehicle. 
The  mission  profile  of  TOPEX  is  to  study  oceanic  sea  surface  topography  using 
instruments  including  an  advanced  radar  altimeter,  a  multi- frequency  radiometer, 
a  laser  retro-reflector  array,  a  TRANET  beacon,  and  an  experimental  high-precision 
radiometric  tracking  device.  The  French  will  supply  a  precision  tracking  system  and 
an  additional  solid  state  altimeter  that  will  be  operational  only  5%  of  the  time  to  avoid 
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interference  with  the  primary  altimeter  [6,28,89,69].  Its  proposed  mission  length  is 
three  years  with  orbital  elements  as  stated  in  Table  1.10.  The  semi-major  axis  for 


Table  1.10.  Orbital  Elements  for  TOPEX 


Keplerian  Element 

TOPEX  Value 

semi-major  axis,  a 

7713.3869  km 

mean  eccentricity,  e 

0.0009825 

inclination,  i 

64.606° 

argument  of  perigee,  u 

270.0° 

longitude  of  the  ascending  node,  fl 

139.552° 

TOPEX  is  significantly  higher  than  the  satellites  reviewed  previously.  This  value  was 
chosen  to  decrease  the  amount  of  atmospheric  drag  on  the  satellite  [69].  To  minimize 
altitude  variations  throughout  the  orbit  and  to  maintain  an  eccentricity  less  than 
0.001,  a  frozen  orbit  is  planned.  In  addition,  the  TOPEX  system  will  have  a  repeat 
ground  track  with  a  repeat  interval  of  10  sidereal  days,  or  127  orbits,  which  will  be 
maintained  within  ±1  km  [6,28]. 

The  inclination  is  constrained  to  meet  several  mission  requirements.  The  intersec¬ 
tion  of  ascending  and  descending  ground  tracks  must  be  nearly  orthogonal  to  deter¬ 
mine  two  orthogonal  components  of  the  surface  with  comparable  accuracy.  Parallel 
tracks  could  accomplish  this  but  would  demand  a  repeat  cycle  longer  than  required 
by  the  mission  specifications.  The  TOPEX  orbit  also  needs  a  nodal  precession  rate 
far  from  sun  synchronous  and  not  near  any  integer  nodal  rate  so  mean  surface  topog¬ 
raphy  can  be  distinguished  from  major  tidal  components.  This  constraint  restricts 
the  TOPEX  inclination  to  be  between  62°  and  65°  [32].  By  examining  repeat  ground 
track  orbits  in  this  inclination  region,  the  TOPEX  orbital  elements  were  chosen  as 
stated  above  in  Table  1.10. 
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1.2.8  Earth  Observing  System  (EOS) 


EOS  consists  of  three  integrated  functions:  the  Scientific  Research  Program,  the 
Data  and  Information  System,  and  the  EOS  Space  Measurement  System  (EOSSMS). 
The  EOSSMS  includes  the  EOS  satellite  system  [29].  Conceptual  studies,  begun  in 
1982,  focused  on  an  optimal  four  or  five  satellite  constellation,  but  later  changed  to 
two  series  of  polar  orbit  platforms,  labelled  A  and  B  [19].  These  satellites,  EOS  A 
B,  will  be  treated  as  single  satellite  systems  and  will  not  be  dependent  on  each 
other  for  their  individual  missions.  Each  platform  series  will  consist  of  three  satellites 
with  five  year  lifetimes,  although  proposals  to  separate  the  B  series  into  several  small 
satellites  are  being  discussed  [25].  These  satellites,  put  up  consecutively,  will  give  each 
series  a  fifteen  year  total  mission  length.  The  A  series  will  focus  on  potential  global 
warming  and  other  aspects  of  global  change  while  the  B  series  will  extend  observations 
made  by  the  Upper  Atmospheric  Research  Satellite  (UARS)  and  TOPEX.  The  delay 
in  determining  the  exact  instrument  complement  for  each  series  allows  for  continuing 
research  in  mission  planning  and  available  technology  [72].  The  first  A  satellite  will 
be  launched  in  1998  followed  by  the  first  B  launch  two  and  a  half  years  later.  The 
current  plan  is  to  launch  all  of  the  satellites  from  the  Western  Test  Range  on  Titan 
IV  rockets  [30]. 

Since  the  mission  constraints  call  for  global  coverage  every  one  to  three  days,  a  sun 
synchronous  orbit  was  chosen  with  a  quasi-two  day  repeat  cycle.  Any  difficulty  with 
separating  the  tides  from  mean  surface  topography  associated  with  a  sun  synchronous 
orbit  is  assumed  to  be  removed  by  sea  geoid  measurements  from  TOPEX  [30,19].  The 
orbital  elements  chosen  to  the  present  are  listed  in  Table  1.11.  Observe  that  these 
orbita'  elements  are  very  similar  to  the  LANDSAT  4/5/6  orbital  elements.  The  semi- 
major  axis  for  EOS  is  only  5  km  higher  than  the  LANDSAT  semi-major  axis.  This 
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Table  1.11.  Orbital  Elements  for  EOS  A/B 


Keplerian  Element 

EOS  A/B  Value 

semi-major  axis,  a 

7083.135  km 

inclination,  i 

98.2° 

argument  of  perigee,  u 

90°  or  270° 

equatorial  crossing  time 

1030  hr  local  time 

semi-major  axis  was  chosen  to  meet  both  the  wide  swath-width  and  high  resolution 
instrument  specifications.  The  inclination  was  chosen  to  insure  a  sun  synchronous 
orbit  at  the  given  altitude  using  equation  1.5.  A  sun  synchronous  orbit  was  chosen 
to  avoid  biasing  the  measurements  of  EOS  with  diurnal  and  seasonal  effects  of  the 
sun  [29].  The  crossing  time  was  altered  from  a  1:30  pm  time  that  complemented 
LANDSAT’s  morning  crossing  times  to  a  10:30  am  time  that  improved  observation 
conditions  in  the  humid  regions  of  the  Earth  during  their  growing  season.  This  new 
choice  of  crossing  time  fails  to  observe  peak  surface  heating,  but  the  disadvantage  of 
this  data  lost  is  still  undetermined  [31]. 

The  A  and  B  series  will  be  supplemented  by  future  National  Oceanic  and  Atmo¬ 
spheric  Administration  (NOAA),  European,  and  Japanese  polar  platforms  planned 
in  conjunction  with  the  EOS  satellites.  EOS  will  also  launch  a  third  series  of  polar 
platforms  dedicated  to  a  synthetic  aperture  radar  (EOS  SAR)  as  on  SEASAT.  This 
instrument,  because  of  its  unique  requirements,  could  not  be  included  on  either  the 
A  or  B  platform  series,  and  will  be  launched  in  1999.  This  series  will  be  composed  of 
three  individual  five  year  lifetime  satellites  that  will  be  placed  in  orbit  consecutively 
for  a  fifteen  year  total  mission  length.  The  EOS  SAR  will  be  in  a  slightly  different 
orbit  than  the  A  or  B  series  with  a  lower  semi-major  axis  of  6998.135  km  but  will 
still  retain  the  1:30  pm  equatorial  crossing  time  [30,29]. 
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1.3  Maneuver  Planning  for  Multiple  Satellites 


Maneuver  planning  for  satellite  constellations  is  no  more  than  a  general  ap¬ 
plication  of  single  satellite  theory  to  many  satellites  with  added  constraints.  These 
constraints  might  be  needed  to  ensure  Earth  coverage,  inter-satellite  communication 
limits,  constellation  geometry,  etc.  There  axe  planar  constraints,  ie.  between  satel¬ 
lites  in  the  same  plane.  The  different  velocities  of  the  satellites  at  separate  parts  of 
the  same  orbit  need  to  be  considered  to  preserve  a  certain  maximum  or  minimum 
angle  between  satellites  in  the  same  plane.  There  are  also  inter-planax  c.  v  saints, 
ie.  between  satellites  of  different  planes,  which  are  more  difficult  to  maintain  since 
the  perturbations  in  the  separate  orbits  can  be  varied.  The  inter-planar  constraints 
may  consist  of  maintaining  phasing  between  satellites  in  different  planes.  A  future 
multiple  satellite  system  that  needs  to  address  these  difficulties  is  the  Iridium  project 
that  has  a  nominal  constellation  of  77  satellites.  It  is  obvious  to  the  casual  observer 
that  the  use  of  sun  synchronous,  repeat  ground  track,  and  frozen  orbits  can  be  applied 
to  more  than  Earth  Observation  satellite  constellations. 

These  same  problems  must  be  addressed  in  US  Defense-oriented  multiple  satel¬ 
lite  constellation  such  as  Brilliant  Pebbles.  The  Soviet  Military  EORSAT  system 
is  thought  to  utilize  an  ocean  surveillance  satellite  constellation  of  6  satellites  in  2 
planes  to  detect,  identify,  and  track  US  and  Allied  naval  forces  [50].  The  two  planes 
are  separated  by  172°  at  the  equator.  The  satellites  have  a  mean  altitude  near  420  km 
and  their  orbits  are  inclined  to  65°.  All  of  the  satellites  are  phased  to  follow  the  same 
repeat  ground  track  cycle  that  lasts  three  days,  or  46  orbit  revolutions.  This  Soviet 
constellation  is  an  example  of  how  the  maintenance  of  a  repeat  ground  track  can 
maintain  the  phasing  needed  by  a  constellation  system.  The  Soviet  EORSAT  system 
lias  had  difficulties  with  keeping  the  correct  number  of  satellites  up  and  keeping  the 
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satellites  in  the  repeat  ground  track  because  of  satellite  errors  in  the  maintenance 
burns.  This  lack  of  accuracy  in  maintenance  burns  lessens  the  synchronicity  of  the 
constellation  [50]. 

1.3.1  Iridium 

The  Iridium  project  sponsored  by  Motorola  Satellite  Communications  is  planned  to 
be  a  commercial,  cellular  communications  system  that  is  satellite  based  for  worldwide 
coverage.  This  system  is  designed  to  complement  the  existing  terrestrial  cellular 
telephone  system.  The  constellation  composed  of  77,  small,  320  kg,  smart  satellites 
will  be  internetted  to  ‘cover’  the  Earth  with  cells  to  allow  channels  to  be  reused  many 
times.  The  portable  units  will  communicate  with  the  constellation  within  a  direct 
line  of  site  on  or  above  the  Earth’s  surface. 

These  77  satellites  will  reside  in  seven  separate  polar  orbital  planes  to  provide 
full  Earth  coverage  with  a  minimum  number  of  satellites  determined  by  Adams  and 
Rider  [1].  These  orbits  are  not  sun  synchronous,  but  to  minimize  orbit  maintenance  a 
frozen  orbit  along  with  a  repeat  ground  track  requirement  are  prescribed,  but  the  final 
decision  is  still  to  be  determined  (TBD).  As  in  the  LANDSAT  4  and  5  projects,  the 
satellites  of  a  constellation  may  be  treated  individually  by  maintaining  each  satellite 
in  a  repeat  ground  track  and  frozen  orbit.  Orbit  inclination  maintenance  will  preserve 
the  ground  track  at  higher  latitudes  and  will  maintain  the  nodal  period  needed  for  a 
repeat  ground  track.  Semi-major  axis  maintenance,  within  some  prescribed  bound¬ 
aries,  will  preserve  the  orbit  that  each  satellite  is  expected  to  follow.  In  addition,  the 
phasing  between  the  satellites  in  the  same  plane  and  other  planes  will  also  be  upheld 
with  maintenance  of  each  individual  repeat  ground  track.  The  frozen  orbit  concept 
would  allow  for  passive  maintenance  of  the  eccentricity  and  argument  of  perigee  of 
each  individual  satellite.  If  all  the  satellites  follow  this  method,  the  interplanar  con- 
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straints  may  be  easier  to  uphold  since  orbit  maintenance  would  only  have  to  focus  on 
the  ground  track  and  not  how  the  other  planes  are  perturbed.  Thus  by  keeping  the 
individual  ground  tracks  intact  for  each  orbital  plane,  inter-planar  angles  would  be 
preserved  without  any  complicated  comparison  between  the  planes.  This  may  not  be 
optimal  in  terms  of  stationkeeping  maneuvers  but  this  possible  solution  does  require 
some  examination. 

The  Iridium  satellites  in  the  odd  numbered  planes  are  in  phase  with  each  other 
and  half  out  of  phase  with  those  in  the  even  numbered  planes.  Most  of  the  planes 
are  co-rotating,  ie.  moving  in  the  same  direction  about  the  Earth.  Only  the  seam 
between  the  first  and  seventh  planes  is  counter-rotating.  The  co-rotating  planes  will 
be  separated  by  27°  and  the  counter-rotating  planes  separated  by  17.5°.  The  distance 
is  smaller  between  the  first  and  seventh  planes  because  Earth  coverage  is  more  difficult 
to  achieve  with  counter-rotating  planes.  These  satellites  will  all  be  in  polar  orbits 
with  a  mean  semi-major  axis  of  7143  km  so  that  the  local  elevation  of  a  satellite  to  a 
portable  user  will  be  greater  than  10°. 

The  maintenance  of  this  satellite  constellation  will  be  difficult,  as  discussed  above, 
because  of  the  many  constraints  in  and  between  the  planes  in  addition  to  the  cellular 
set  up  of  the  communication  linkage  [56]. 

1.4  Summary 

SEASAT,  GEOSAT,  and  TOPEX  are  repeat  ground  track  and  frozen  orbit  missions 
that  do  not  utilize  the  sun  synchronous  condition.  Conversely,  EOS,  LANDSAT  4/5, 
NROSS  and  the  European  Space  Agency  (ESA)  ERS-1  all  utilize  the  three  concepts 
mentioned  above.  This  is  summarized  in  Table  1.12. 

This  thesis  will  focus  on  maneuver  planning  for  a  sun  synchronous,  frozen  orbit, 
repeat  ground  track  orbit  applied  to  the  LANDSAT  6  mission  constraints.  Instead 
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Table  1.12.  Satellite  Orbit  Summary 


Satellite 

System 

Repeat 

Ground  Track 

Sun 

Synchronous 

Frozen 

Orbit 

SEASAT 

YES 

NO 

YES 

LANDSAT  4/5 

YES 

YES 

YES 

NROSS 

YES 

YES 

YES 

GEOS  AT  ERM 

YES 

NO 

YES 

ERS-1 

YES 

YES 

YES 

TOPEX 

YES 

NO 

YES 

EOS 

YES 

YES 

YES 

Iridium 

TBD 

NO 

TBD 

of  using  more  common  propagation  methods  as  in  past  satellite  systems,  this  thesis 
takes  a  new  approach  to  maneuver  planning  by  suggesting  the  employment  of  a  high 
precision  mean  element  orbit  propagator,  called  the  Draper  Semianalytic  Satellite 
Theory  (DSST)  to  predict  how  the  orbit  changes  from  the  nominal  elements  over 
time. 

This  propagation  method  is  not  only  accurate  and  computationally  efficient,  but 
it  is  also  singularity  free  through  the  usage  of  equinoctial  orbital  elements  rather  than 
the  more  common  Keplerian  elements.  The  DSST  Standalone  Orbit  Propagator  is 
small  and  its  modern  architecture  allows  it  to  be  portable  to  many  computer  systems 
including  IBM  and  CDC  Mainframes,  a  VAX,  an  IBM  PC,  a  Sun  SPARCstation,  and 
a  Macintosh  PC. 

Chapter  2  reviews  the  mean  element  theory  and  program  flow  used  in  the  DSST 
Orbit  Propagator.  Chapter  3  discusses  an  existing  software  tool,  MEANELT,  that 
applies  the  DSST  propagation  techniques  to  maneuver  planning  for  maintaining  the 
nominal  orbit.  Chapter  4  states  the  specific  algorithms  necessary  for  sun  synchronous, 
repeat  ground  track,  and  frozen  orbit  maneuver  planning  and  specifies  the  impulsive 
and  finite  burn  models  used  in  the  LANDSAT  6  mission  software.  Chapter  5  discusses 
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ground  track  motion  forecasting  using  solar  flux  and  geomagnetic  activity  predictions 
and  the  confidence  interval  of  the  density  calculated  form  these  values.  Chapter  6 
asserts  the  results  and  conclusions  of  this  thesis  and  gives  recommendations  for  fu¬ 
ture  work.  A  review  of  the  Keplerian  and  equinoctial  elements  sets  can  be  found  in 
Appendix  A.  Appendix  B  discusses  the  porting  of  the  Draper  Semianalytic  Satellite 
Theory  Standalone  Propagator  from  the  IBM  mainframe  to  other  computing  envi¬ 
ronments,  ie.  VAX,  Sun  SPARCstation,  IBM  PC,  and  Macintosh  PC. 
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Chapter  2 

Review  of  Semianalytic  Theory 


To  determine  a  satellite  orbit  precisely  in  time  and  to  insure  that  it  meets  the 
mission  specifications,  the  orbit  must  be  propagated  through  time  using  perturbation 
methods  from  two  basic  categories:  Special  and  General  Perturbation  Theories.  The 
first,  called  the  Numerical  Method  or  the  Special  Perturbation  Theory,  constructs  the 
perturbed  orbit  through  direct  integration  of  the  equations  of  motion.  This  method 
is  very  precise  but  computationally  inefficient;  to  retain  accuracy,  the  high  frequency 
perturbations  constrain  the  integration  step  size  to  be  very  small.  In  addition,  the 
large  number  of  steps  taken  increase  the  truncation  and  round-off  error  in  many 
orbit  determination  applications.  The  Special  Perturbation  Method  does  not  provide 
general  physical  insight  into  orbit  dynamics.  A  Cowell  propagator  is  an  example  of  a 
Special  Perturbation  Theory.  On  the  other  hand,  the  General  Perturbation  Theory 
or  the  Analytic  Model,  transforms  the  equations  of  motion  into  exact  differentials 
that  can  be  analytically  integrated.  To  realize  these  closed  form  analytic  expressions, 
the  Keplerian  orbit  including  all  perturbations  is  generalized  through  assumptions,  ie. 
simplified  perturbation  models,  approximations,  ie.  small  terms  in  series  expansions 
are  neglected,  and  restricted  ranges  of  theoretical  validity.  This  method  is  not  very 
accurate  because  of  the  many  simplifications,  but  it  is  computationally  efficient  after 
the  algorithms  have  been  defined.  If  any  additional  perturbations  are  desired,  another 


tedious  derivation  must  take  place  to  construct  the  new  analytic  expressions  [4,23]. 

An  alternative  to,  or  combination  of,  these  two  classes  of  perturbation  methods 
is  a  semianalytic  theory  that  is  both  accurate  and  computationally  efficient.  Poten¬ 
tial,  or  conservative,  perturbations  are  put  in  Lagrangian  Variation  of  Parameters 
(VOP)  form  and  non-potential,  or  nonconservative,  perturbations  are  put  in  Gaus¬ 
sian  Variation  of  Parameters  form.  The  long  period  and  secular  components  of  the 
perturbations  are  separated  from  the  short  period  components  of  satellite  motion.  In 
the  mean  element  theory  used  in  this  thesis,  this  is  done  through  the  application  of 
the  Generalized  Method  of  Averaging  to  the  VOP  equations  of  motion.  The  simple 
conservative  perturbations  are  analytically  averaged.  These  conservative  perturba¬ 
tions  could  be  numerically  averaged,  but  this  method  is  not  computationally  efficient, 
thus  it  is  reserved  for  the  more  complex,  non-conservative  perturbations.  Through 
these  averaged  equations  of  motion,  integration  of  the  mean  element  (long  period  and 
secular)  motion  is  achieved  efficiently  by  applying  large  step  sizes.  The  magnitude  of 
this  step  size  is  constrained  by  the  magnitude  of  the  next  higher  frequency  oscilla¬ 
tions  [87,23].  In  this  thesis,  the  step  size  is  one  day  or  larger  [27].  Accuracy  is  achieved 
by  adding  in  the  short  periodic  terms  at  the  output  time.  In  the  Draper  Semianalytic 
Satellite  Theory  Standalone,  these  short  periodic  terms  are  determined  on  a  grid  and 
if  the  output  time  is  not  on  the  grid,  the  short  periodic  terms  are  interpolated.  The 
only  difficulty  with  the  mean,  semianalytic  theory  is  that  it  requires  a  procedure  to 
determine  the  initial  mean  elements  to  start  the  integration  at  the  same  osculating 
elements  as  a  numerical  integrator.  This  problem  only  arises  when  the  semianalytic 
theory  is  to  be  compared  to  a  numerical  integration  of  the  same  satellite  orbit.  In  nor¬ 
mal  orbit  determination  operations,  convergence  of  the  semianalytic  method  rapidly 
eliminates  errors  in  the  mean  initial  conditions.  For  comparison  testing,  the  initial 
mean  elements  are  determined  through  a  least  squares  fit  of  the  semianalytic  theory 
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to  the  desired  osculating  trajectory.  This  mean  initial  state  could  also  be  achieved 
by  numerically  averaging  the  osculating  state  over  one  orbit  [87,23].  A  more  rigorous 
mathematical  summary  is  stated  below. 

There  is  ongoing  research  into  new  methodologies  for  the  development  of  analytical 
and  semianalytical  theories.  These  new  methods  explore  the  use  of  algebraic  software 
and  hardware  to  decrease  the  work  in  deriving  the  analytic  expressions  [24], 

2.1  Mathematical  Summary 

2.1.1  Variation  of  Parameters 

Newton’s  second  law  may  be  applied  to  the  equations  of  motion  of  two  point  mass 
particles,  ie.  the  Earth  and  the  satellite,  to  determine  the  fundamental  differential 
equation  of  the  two-body  Keplerian  orbit  [4]: 

r  +  ^r  =  0  (2.1) 

where  r  is  the  position  vector  of  the  satellite  with  magnitude  r,  yu  is  the  gravitational 
parameter  of  the  Earth,  and  over  dots  denote  time  differentiations.  The  solution  to 
this  equation  leads  to  six  constants  of  integration: 

C  =  f  C.  c2  c3  c4  c5  Ce  ]T  (2.2) 

which  can  be  represented: 

r  =  r(c,  t)  (2.3) 

r  =  v(c  ,t)  =  ^  (2-4) 

where  r  and  v  represent  the  transformations  between  the  position/ velocity  vectors 
and  the  orbital  elements  of  c.  The  most  common  of  these  orbital  elements  are  the 
Keplerian  orbital  elements,  but  a  variation  of  these  are  the  equinoctial  elements;  both 
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described  in  Appendix  A.  These  orbital  elements  describe  the  two  body,  osculating, 
orbit  that  would  exist  in  a  system  only  affected  by  the  law  of  gravitation  with  the 
bodies  involved  defined  as  point  masses.  But,  a  satellite  orbiting  the  Earth  is  not 
influenced  only  by  the  gravitational  attraction  of  the  Earth,  meaning  that  the  orbit  is 
also  affected  by  perturbations  to  this  two- body  orbit.  This  is  denoted  by  the  following 
equation: 

r  +  =  ad  (2.5) 

where  ad  is  the  perturbing  acceleration  acting  on  the  satellite.  If  the  perturbing 
acceleration,  a^,  is  very  small  compared  to  the  primary  gravitational  acceleration,  the 
solution  to  Equation  2.5  would  closely  approximate  the  solution  to  the  unperturbed 
case  in  Equation  2.1.  Under  the  additional  perturbations,  ad,  the  orbital  elements  in 
c  change  slowly  over  time.  The  purpose  of  the  Method  of  Variation  of  Parameters 
(VOP)  is  to  construct  differential  equations  which  give  the  motion  of  these  orbital 
elements  over  time. 


Lagrangian  Variation  of  Parameters 


The  Lagrange  VOP  equations  express  the  perturbing  acceleration  as  the  gradient 
a  conservative  disturbing  function,  R: 


*d  = 


[<«1 

T 

SR 

SR 

SR  ] 

1 - 

Or> 

*1 

1 _ 

Sr  i 

Sr  2 

S  r3  J 

(2.6) 


where  r i,  r2,  and  r3  are  the  coordinates  of  the  position  vector,  r.  Differencing  Equa¬ 
tion  2.3  and  allowing  the  orbital  elements  of  c  to  vary  with  time  gives: 


v 


<5r  dc 
6c  dt 


(2.7) 
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By  substituting  the  latter  part  of  Equation  2.4  into  the  above  equation,  the  first  three 
constraints  on  the  orbital  elements  c  are  derived: 


6r  dc 
—  —  =  0 
oc  dt 


Next  by  differencing  Equation  2.4,  an  expression  for  the  second  derivative  of  the 
position  vector,  r,  is  determined  and  can  be  substituted  into  the  perturbed  equation 


of  motion  giving: 


<52r  6v  dc  n  6R 

6t2  +  6c  dt  +  t3*  6 r 


Since  c  is  the  solution  for  the  unperturbed  case,  Equation  2.1  is  restated  in  the 
following  form  and  substituted  into  Equation  2.9  to  give  the  final  three  conditions  on 


the  orbital  elements  c: 


<52r  /x 

+  -^r  =  0 


6v  dc 
6c  dt 


(2.10) 

(2.11) 


The  six  constraints,  or  differential  relations,  of  the  orbital  elements,  c,  of  Equations  2.8 


and  2. 1 1  can  be  placed  in  matrix  form: 


ft  = 


By  pre-multiplying  this  equation  by  the  6x6  matrix: 


and  using  the  chain  rule  of  partial  differentiation: 

6R  _  6R6r 
6c  6r  6c 


(2.12) 


(2.13) 


(2.14) 


a  convenient  and  more  familiar  expression  for  the  transformation  between  the  vari¬ 


ables  r,  r  and  the  orbital  elements,  c  is  provided: 


dc 

L  ~~  — 
dt  6c 


(2.15) 


where  the  skew-symmetric  Lagrangian  matrix  is  defined: 
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(2.16) 

(2.17) 


To  achieve  the  orbital  element  rates  directly  from  Equation  2.15,  the  inverse  of  the 
Lagrangian  matrix  L  must  be  determined.  The  proof  given  by  de  Lafontaine  [23] 
shows  that  the  skew  symmetric  Poisson  matrix,  P,  is  the  negative  inverse  of  L.  This 
gives  the  Lagrangian  VOP  equations  of  motion  [4]: 

dirpTt  ^ 


where  the  Poisson  bracket,  P  is  defined: 
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Gaussian  Variation  of  Parameters 


(2.19) 


The  Gaussian  form  of  Variation  of  Parameters  is  developed  with  no  assumptions 
made  concerning  the  perturbing  accelerations.  This  differs  from  the  Lagrangian  form 
where  the  perturbing  acceleration  is  defined  as  the  gradient  of  some  disturbing  func¬ 
tion,  R.  The  Gaussian  VOP  equations  of  motion  are  used  for  perturbations  that  can 
not  be  expressed  in  terms  of  R ,  ie.  non-conservative  forces  like  atmospheric  drag.  By 
no  longer  allowing  the  substitution  of  Equation  2.6,  Equation  2.18  is  now  stated: 

i  t 


dc  _  pT 

dt 


6  r 
8c 
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(2.20) 


By  expanding  the  Poisson  bracket,  the  following  equation  is  achieved: 
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(2.21) 
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Since  the  components  of  the  position  and  velocity  vectors  are  considered  independent, 
the  following  equations  hold  true: 


8  r 

<5r 

8r 

8v 


=  I 
=  O 


(2.22) 

(2.23) 


where  I  is  the  identity  matrix  and  0  is  the  zero  matrix.  These  equations  may  be  sub¬ 
stituted  into  Equation  2.21  to  result  in  the  final  form  of  the  Gaussian  VOP  equations 
of  motion  [4]: 


dc  8c 
dt  8\Ad 


(2.24) 


This  Gaussian  form  may  be  applied  to  determine  the  rates  of  change  of  the  sin¬ 
gularity  free  equinoctial  element  set  [57]: 
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(2.26) 

(2.27) 

(2.28) 

(2.29) 

(2.30) 


The  parameter,  G,  is  determined  by  the  following  equation: 


G  =  na2V  1  —  h2  —  k2 


The  satellite  mean  motion,  n,  is  calculated  with  the  equation: 


n  =  r  — 


(2.31) 


(2.32) 
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The  variables,  AT i ,  Y\ ,  X\,  and  Yy,  are  the  position  and  velocity  coordinates  of  the 
satellite  in  the  equinoctial  orbit  frame.  They  are  determined  with  the  following  equa¬ 
tions: 


=  a[(l  —  h2/3)  cos  F  +  hk(3  sin  F  —  k] 
=  a[(l  —  k2/3)sinF  +  hk/3  cos  F  —  h] 

77 

—  - [hk/3  cos  F  —  (1  —  h2j3)  sinF] 


(2.33) 

(2.34) 

(2.35) 


-  [( 1  —  k2(3 )  cos  F  —  hk/3  sin  F] 


(2.36) 


where  the  parameter,  /?,  is  equated: 


1  +  v/1  -  h2  -  k2 


(2.37) 


The  eccentric  longitude,  F,  can  be  determined  through  an  iterative  solution  of  the 


following  equation: 


A  =  F  +  h  cos  F  —  k  sin  F 


(2.38) 


The  unit  vectors  /  and  g  both  lie  in  the  orbit  plane  while  the  unit  vector  w  is 
perpendicular  to  the  orbit  plane.  This  orthogonal  coordinate  system  is  determined 


through  the  transformation: 

2  1  —  p2  +  q2  2  pql  2 p 

1  f  9  w]  =  — — 2-7-2  2 pq  /(I  +  p2  -  q2)  -2 q 

P  q  —2  Ip  2  q  I(l—p2  —  q2) 


(2.39) 


The  position  and  velocity  vectors  are  now  determined  in  this  coordinate  system: 


r  =  XJ  I  Yxg 
v  -  XJ  \  Yxg 


(2.40) 

(2.41) 


Since  the  DSST  Standalone  uses  both  the  Lagrangian  and  the  Gaussian  VOP 
methods,  the  perturbing  acceleration  is  separated  into  conservative  and  non-conservative 
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parts.  This  loads  to  a  general  form  of  the  VOP  equations  of  motion  [23]. 

.  .  T 

l  6t 

&d  **d.\c  "f"  *^dc  *^d-NC  "f- 


<5r 


dc  6  R  6c 

di  =  P  foadNC 


(2.42) 

(2.43) 


2.1.2  Generalized  Method  of  Averaging 

The  equations  of  motion  in  both  the  Lagrangian  and  Gaussian  VOP  forms  include 
the  secular  and  the  long  and  short  periodic  perturbations  that  disturb  the  orbit 
from  a  two-body  Kepler ian  orbit.  The  Generalized  Method  of  Averaging  eliminates 
the  short  periodic  terms  from  these  expressions  to  form  mean  element  equations  of 
motion.  These  mean  element  equations  may  be  integrated  using  a  larger  time  step, 
thus  decreasing  the  amount  of  computational  time.  The  short  periodic  terms  can  be 
added  in  at  the  requested  time  to  recall  the  accurate  osculating  state  of  the  satellite. 

To  determine  the  mean  element  equations  of  motion,  the  orbital  elements  must 
be  separated  into  a  group  of  slow  osculating  orbital  elements,  c,  and  a  fast  variable, 
c6.  This  separation  in  the  VOP  equations  of  motion  can  be  shown  functionally  in  the 
following  form  [60]: 

^  =  eG(c,  cq)  (2.44) 

“  =  /i(c)  +  eG6(c,ce)  (2.45) 


The  functions  G  and  G6  represent  the  perturbing  accelerations  from  the  right  hand 
side  of  the  equations  of  motion  while  the  function  h  represents  the  Keplerian  motion  of 
the  fast  parameter,  see  Equation  2.43.  A  near  identity  first  order  transformation  from 
the  mean  to  osculating  elements  is  determined  with  the  following  expressions  [88]: 


c  =  c  I  e7?(c,c6) 
C(j  Co  i  £?76(d,co) 
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(2.46) 

(2.47) 


where  the  over  bar  denotes  the  averaged  elements  and  the  functions  er)  and  erje  denote 
the  short  periodic  functions,  and  are  the  first  order  terms  in  an  asymptotic  series  for 
the  original  dynamics.  The  rates  of  the  mean  elements  can  not  depend  on  the  variables 
e6  or  c6  because  these  are  fast  changing  values  and  are  eliminated  in  the  averaging. 
Therefore,  the  mean  element  rates  can  be  written: 

c  =  eA(c)  (2.48) 

he  =  h(c)  4-  e,46(c)  (2.49) 

where  again  the  variable  e  denotes  the  small  magnitudes  of  the  disturbing  accelera¬ 
tions,  A  and  A&.  The  mean  element  rates  and  short  periodic  functions  are  determined 
by  first  assuming  an  asymptotic  matching  between  the  original  VOP  equations  and 
equations  resulting  from  the  differentiation  of  Equations  2.46  and  2.47.  This  match¬ 
ing  is  usually  done  to  first  and  second  order  in  the  variable  e.  The  expressions  for  the 
mean  element  rates  and  short  periodic  functions  are  realized  with  the  the  substitution 
of  the  mean  element  rates,  Equations  2.48  and  2.49,  into  the  asymptotic  expansion. 
One  important  constraint  is  that  the  short  periodic  functions,  r/  and  %,  are  periodic 
in  the  fast  variable,  Cg,  with  a  period  of  2-7t;  ie.  the  following  constraints  apply  [88]: 

€77(c,  c6  +  27t)  =  er)(c,  eg)  (2.50) 

f  trj(c,c6)dce  =  0  (2-51) 

Jo 

This  assumption  allows  simplification  of  the  mean  element  rates  so  that  first  order 
coefficients,  A,  are  calculated  [60]: 

A,  =  ^  Gi(c,c6)dc6  (2.52) 

where  the  integer  i  is  from  1  to  6. 
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2.1.3  Mean  Element  Equations  of  Motion 

For  conservative  forces,  like  the  Earth’s  geopotential  zonal  harmonics,  the  mean  ele¬ 
ment  equations  of  motion  are  found  through  an  analytical  averaging  of  the  Lagrangian 
VOP  equations  of  motion  over  the  fast  variable  for  one  full  orbital  period  [60,88]: 

1  /-2w  f)R 

eA'  =  2^i  P  JZdc«  (2'53) 

But,  the  Poisson  matrix,  P,  depends  only  on  the  slowly  varying  elements,  c,  so  it 
may  be  pulled  through  the  integral. 

1  f  2ir  f.  £? 

tAi  =  PTi~  /  -rdct  (2.54) 

27t  Jo  oc 

By  interchanging  the  order  of  the  integral  and  the  differential,  the  right  hand  sides  of 
the  mean  element  equations  of  motion  are  expressed  in  terms  of  the  averaged  potential 
and  the  partial  derivatives  of  the  averaged  potential: 

tA‘  =  pTrc[hfmc*}  (255) 

The  average  of  the  potential  and  the  derivatives  of  the  averaged  potential  have  been 
constructed  for  the  zonal  harmonics,  tesseral  harmonics,  and  the  third  body  point 
mass  effects.  Detailed  mathematical  expressions  for  the  averaged  geopotential  zonal 
harmonics  can  be  found  in  reference  [16],  the  effects  of  Jf  *n  reference  [8],  the  averaged 
gcopotential  tesseral  harmonics  (resonance)  can  be  found  in  references  [76,14],  and 
the  averaged  effects  of  solar-lunar  point  masses  in  reference  [16]. 

For  non-conservative  perturbations,  like  atmospheric  drag  and  solar  radiation,  the 
Draper  Scmianalytic  Theory  applies  the  general  first-order  averaging  method  to  the 
singularity  free  equinoctial  element  equations  of  motion  described  in  Equations  2.25 
through  2.30.  The  Method  of  Averaging  is  still  applied  but  now  the  Gaussian  VOP 
equations  of  motion  are  averaged  using  a  Gaussian  quadrature  method  [60,88]. 
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In  the  Keplerian  element  set,  the  fast  element  is  represented  by  the  mean  anomaly 
and  in  the  equinoctial  element  set  it  is  the  mean  longitude.  To  uniformize  the  inte¬ 
grand  and  increase  computational  efficiency,  the  integration  variable  is  transformed 
from  the  mean  longitude,  A,  to  the  eccentric  longitude,  F,  with  the  following  expres¬ 
sion  [57]: 

dp 

— =-  =  (1  —  kcosF  —  /isinF)-1  (2.56) 

dX 

2.1.4  Short  Periodic  Perturbations 

If  only  the  mean  element  equations  of  motion  are  integrated,  accuracy  is  lost. 
Therefore,  the  short  periodic  terms  must  be  recalled  and  added  in  to  determine  the 
true  osculating  state  at  the  requested  output  time.  If  the  short  periodic  terms  are 
calculated  at  every  output  time,  there  is  a  loss  of  computational  efficiency  if  the 
desired  output  times  are  at  many  closely  spaced  points  in  time.  The  Draper  Semian- 
alytical  Satellite  Theory  Standalone  (DSST  Standalone)  determines  the  coefficients 
of  the  short  periodic  Fourier  series  expansions  on  a  grid.  This  grid  is  closely  tied  to 
the  integration  of  the  mean  element  equations  of  motion.  If  the  request  time  is  in 
between  the  grid  points,  the  short  periodic  coefficients  are  determined  with  a  three 
point  Lagrangian  interpolator  formula.  If  the  request  time  is  not  spanned  by  coeffi¬ 
cients  already  calculated,  the  program  determines  the  correct  interpolation  interval 
and  calculates  the  coefficients  and  uses  these  to  interpolate  the  short  periodics  at  the 
request  time.  This  formulation  allows  a  high  execution  speed  for  many  closely  spaced 
request  times  [27]. 

As  stated  above,  the  short  periodics  are  determined  through  a  Fourier  series  ex¬ 
pansion  in  the  fast  variable  with  C,  and  D,  representing  the  series  coefficients  [88] : 

N 

tr](c,ce)  =  E‘C.(c)  sin  icQ  —  eDj(c)  cos  ic$  (2-57) 

t=i 
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(2.58) 


N 

er)6(c ,  ce)  =  e(^6,i(c)  sin*C6  -  e-D6,»( c)  cos  ic^ 

i- 1 

In  the  DSST  Standalone,  the  fast  variable,  eg  is  either  the  mean  longitude,  A,  eccen¬ 
tric  longitude,  F,  or  true  longitude,  L  [27].  The  choice  of  which  longitude  to  use  as 
the  expansion  variable  is  different  for  various  perturbations.  The  central  body  zonal 
harmonics  are  in  a  closed  form  expression  if  the  expansion  variable  is  the  true  lon¬ 
gitude,  L.  The  same  is  true  for  the  slow  moving  third  body  perturbations  with  the 
eccentric  longitude,  F.  A  more  detailed  examination  of  the  mathematical  expressions 
is  given  for  the  zonal  harmonics  in  references  [51,81],  third  body  short  periodic  effects 
in  references  [80,51],  and  the  tesseral  harmonics  in  references  [51,77,17,76,14]. 

The  central  body  gravitational  sectoral  and  tesseral  harmonics  short  periodics  are 
determined  through  two  expansions.  The  first  expansion  is  in  the  longitude  of  the 
central  body’s  prime  meridian,  9,  which  is  measured  from  the  equinoctial  origin  of 
longitudes  in  the  satellite  orbit.  These  periodics  require  a  separate  expansion  since 
they  are  of  medium  period,  ie.  T  times  the  rotational  period  of  the  central  body 
where  m  is  a  small  integer.  The  second  expansion  is  in  both  the  mean  longitude,  A, 
and  the  longitude  of  the  prime  meridian,  9.  This  expansion  leads  to  short  periodic 
terms  with  a  much  higher  frequency  than  the  satellite  period.  This  double  expansion 
may  result  in  some  very  long  period  erms  resulting  from  tesseral  resonance.  These 
long  periodic  terms  are  removed  from  the  short  periodic  expansion  and  added  to  the 
averaged  equations  of  motion. 

For  third  body  point  mass  perturbations,  there  is  a  choice  of  two  formulations  to 
determine  the  associated  short  periodics.  The  first  formulation  employs  an  expansion 
in  the  eccentric  longitude  of  the  satellite,  F  [80].  The  second  formulation  employs 
both  the  mean  longitude  of  the  satellite,  A,  and  the  mean  longitude  of  the  third  body, 
A',  as  trigonometric  variables.  This  last  formulation  is  particularly  useful  for  high 


57 


altitude  satellites  whose  mean  motion  is  not  large  relative  to  the  mean  motion  of  the 
perturbing  third  body.  An  Earth  Observation  satellite  has  a  relatively  short  period 
and  the  third  body  can  be  assumed  to  have  a  constant  position  over  the  satellite 
orbital  period  [27]. 

2.2  Draper  Semianalytic  Satellite  Theory  Standalone 
Orbit  Propagator  (DSST  Standalone) 

The  Draper  Semianalytic  Satellite  Theory  Standalone  Orbit  Propagator  is  a 
descendant  of  the  Goddard  Trajectory  Determination  System  (GTDS)  Averaged  Or¬ 
bit  Generator.  The  GTDS  program  was  modified  by  Draper  Laboratory  to  include 
more  perturbations,  interpolator  choices,  models  for  the  short  periodic  variations,  and 
partial  derivatives  of  the  current  mean  elements  with  respect  to  the  solve  for  param¬ 
eters  [27],  The  DSST  Standalone  is  a  straight  orbit  propagation  program  version  of 
the  Draper  GTDS  program  which  is  a  full  orbit  determination  system.  The  DSST 
Standalone  is  portable  to  different  host  programs  and  to  different  computing  envi¬ 
ronments.  This  thesis  will  focus  on  the  DSST  Standalone  because  of  its  application 
to  a  maneuver  planning  simulation. 

The  Draper  Semianalytic  Satellite  Theory  Standalone  Orbit  Propagator  initializes 
and  propagates  orbital  elements,  or  ephemeris  data,  over  a  period  that  is  maneuver 
free.  If  a  maneuver  takes  place,  the  DSST  Propagator  reinitializes  the  ephemeris 
data  and  continues  the  propagation  scheme.  The  driver  program,  SEMIANAL, 
utilizes  five  basic  subroutines  shown  in  Figure  2.1.  Subroutine  SETELM  initializes 
the  epoch  time,  epoch  elements,  and  the  duration  of  the  propagation.  INTANL 
initializes  the  mean  clement  propagator.  A  call  to  BEGANL  sets  the  direction  of 
integration  and  starts  the  integrator  at  the  epoch  time.  If  the  user  requests  output 
at  specified  times,  subroutine  ORBANL  will  propagate  to  that  time  and  output 
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SEMIHNHNL 


Figure  2.1.  DSST  Driver  SEMIANAL 


the  ephemeris.  Subroutine  RESANL  restarts  the  integration  after  an  impulsive 
maneuver  takes  place.  If  the  maneuver  is  not  modeled  as  impulsive,  calls  to  INTANL 
and  BEGANL  will  restart  the  integrator  to  set  the  new  epoch  ephemeris  to  the 
orbital  elements  after  the  maneuver.  Calls  to  ORBANL  will  continue  until  the  end 
of  the  propagation  period  is  reached. 

2.2.1  SETELM 

Subroutine  SETELM  sets  the  initialization  parameters  needed  by  the  DSST 
Propagator  to  start  the  integration.  The  parameters  which  should  be  set  in  this 
subroutine  and  their  definitions  are  listed  in  Table  2.1. 

2.2.2  INTANL 

Subroutine  INTANL  is  called  at  the  beginning  of  an  orbit  propagation  to 
initialize  the  integrator  parameters  and  force  models.  The  parameters  input  into  this 
subroutine  are  listed  in  Table  2.2. 

Subroutine  INTANL  calls  three  subroutines  to  initialize  the  integrator  parame¬ 
ters  and  force  models.  This  is  shown  in  Figure  2.2.  Subroutine  SETANL  first  calls 
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Table  2.1.  SETELM  Output 


Variable 

Definition 

ELMINT 

Mean  orbital  elements  at  epoch. 

If  Keplerian,  angular  values  are  in  deg 
a,  e,  i,  12,  u,  M 

If  Equinoctial,  angular  values  are  in  rad 
u,  /i,  Aj,  p,  <y,  A 

RETINT 

Retrograde  factor  if  ELMINT 
is  equinoctial 

ITYPE 

Type  of  input  orbital  elements 
=  1,  Cartesian  (position,  velocity) 

=  2,  Keplerian 
=  3,  Equinoctial 

ICOORD 

Coordinate  system  of  orbital  elements 
=  1,  Mean  1950 
=  2,  True  of  Date 

IOSCU 

Are  elements  osculating  or  mean? 

=  1,  osculating 
=  2,  mean 

YMDINT 

Calendar  date  of  epoch 

YYYYMMDD. 

HMSINT 

Time  of  epoch 

HHMMSS.SSSSS 

YMDEND 

Calendar  date  of  end 

YYYYMMDD. 

HMSEND 

Time  of  end 

HHMMSS.SSSSS 

PRTSTP 

Print  interval  in  UTC  seconds 

Figure  2.2.  Subroutine  INTANL 
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Table  2.2.  INTANL  Input 


Variable 

Definition 

ELMINT 

Mean  orbital  elements  at  epoch. 

If  Keplerian,  angular  values  are  in  deg 
a,  e,  i,  fl,  u,  M 

If  Equinoctial,  angular  values  are  in  rad 
a,  h ,  k ,  p,  q,  A 

RETRO 

Retrograde  factor  if  ELMINT 
is  equinoctial 

ITYPE 

Type  of  input  orbital  elements 
=  1,  Cartesian  (position,  velocity) 

=  2,  Keplerian 
=  3,  Equinoctial 

ICOORD 

Coordinate  system  of  orbital  elements 
=  1,  Mean  1950 
-  2,  True  of  Date 

IOSCU 

Are  elements  osculating  or  mean? 

=  1,  osculating 
=  2,  mean 

YMDINT 

Calendar  date  of  epoch 

YYYYMMDD. 

HMSINT 

Time  of  epoch 

HHMMSS.SSSSS 
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SETFRC  to  initialize  the  force  model  parameters.  Then  a  call  to  SETAVR  sets 
the  parameters  relevant  to  the  mean  element  equations  of  motion.  SETANL  then 
calls  SETSP  to  set  the  parameters  related  to  the  model  of  the  short-periodic  pertur¬ 
bations.  Finally,  SETANL  calls  SETDRV  to  set  flags  related  to  the  computation 
of  partial  derivatives.  Then  INTANL  calls  subroutine  ANLINT  to  determine  the 
parameters  that  are  needed  to  initialize  the  integrator.  ANLCRD  is  presently  an 
empty  subroutine;  in  the  future,  ANLCRD  will  process  card  inputs  to  the  program. 

2.2.3  BEGANL 

Subroutine  BEGANL  is  called  at  the  beginning  of  an  orbit  propagation  to 
start  the  integrator  at  the  desired  epoch.  This  is  accomplished  by  setting  the  inte¬ 
gration  direction  and  the  interpolation  time  step  parameters.  The  propagator  is  able 
to  propagate  forwards  or  backwards  in  time,  but  it  may  not  switch  directions  in  the 
middle  of  a  propagation.  It  then  initializes  the  partial  derivatives  at  epoch.  Subrou¬ 
tine  RESANL  is  called  to  start  the  integrator.  BEGANL  must  be  preceded  by  a 
call  to  INTANL.  One  of  the  needed  input  parameters,  DIRINT,  is  input  through 
BEG  ANL’s  argument  listing  and  the  rest  are  input  through  COMMON  blocks  de¬ 
noted  by  surrounding  slashes.  The  input  parameters  are  described  in  Table  2.3.  The 
output  of  BEGANL  is  all  through  COMMON  blocks  and  is  listed  in  Table  2.4. 

2.2.4  ORBANL 

Subroutine  ORBANL  is  called  to  return  the  satellite  state  corresponding  to 
the  requested  time.  The  subroutines  called  by  ORBANL  can  be  viewed  in  Fig¬ 
ure  2.3.  If  the  requested  time  is  not  within  the  interval  of  calculated  coefficients, 
subroutine  STPANL  is  called  to  propagate  the  mean  orbital  elements,  calculate  the 
mean  element  rates  and  determine  the  short  periodic  interpolator  coefficients  for  the 
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Table  2.3.  BEGANL  Input 


Variable 

Definition 

DIRINT 

Direction  of  integration 
=  1,  forward 
=  -1,  backward 

/ANLFIL/ 

NPRINT 

Standard  output  unit 

/ANLPRM/ 

ELMEPO 

Mean  equinoctial  elements  at  epoch 

NSOLVE 

Number  of  solve  for  parameters 

NSTATE 

Number  of  state  solve  for  parameters 

NDYNAM 

Number  of  dynamic  solve  for  parameters 

/AVRHST/ 

IIAFWID 

Half-width  of  interpolator  interval 

STEP 

Integrator  step  size,  negative  if  backwards 

/SPREAL/ 

PVSTEP 

Interval  between  successive 
interpolator  points 

PVWID 

Half-width  of  interpolator  interval 

SPSTEP 

Interval  between  successive 
interpolator  points 

SPWID 

Half-width  of  interpolator  interval 

Table  2.4.  BEGANL  Output 


Variable 

Definition 

/AVRHST/ 

HAFWID 

Half-width  of  interpolator  interval 

STEP 

Integrator  step  size,  negative  if  backwards 

AVRDRV 

Partials  of  current  mean  elements 
with  respect  to  solve  for  parameters 

/SPREAL/ 

PVSTEP 

Interval  between  successive 
interpolator  points 

PVWID 

Half-width  of  interpolator  interval 

SPSTEP 

Interval  between  successive 
interpolator  points 

SPWTD 

_ i 

Half-width  of  interpolator  interval 
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Figure  2.3.  Subroutine  ORBANL 


new  interval.  Subroutine  ELEMNT  determines  the  short  periodic  terms  that  need 
to  be  added  to  the  mean  elements  to  get  the  osculating  satellite  state  at  the  out¬ 
put  time.  Subroutine  ELEMNT  is  called  much  more  frequently  than  subroutine 
STPANL  since  the  mean  elements  and  the  short  periodics  can  be  interpolated  from 
one  integration  interval’s  coefficients  for  several  output  times. 

If  an  external  event  prematurely  ends  the  propagation,  the  time  of  the  propagation 
stop  is  returned  with  the  satellite  state  and  partial  derivatives  at  that  time.  The  input 
parameters  needed  by  ORBANL  are  listed  in  Table  2.5.  The  satellite  state  at  the 
request  time  is  output  in  part  through  ORBANL’s  argument  listing  and  the  rest 
through  COMMON  blocks.  The  exact  output  is  listed  in  Table  2.6. 

2.2.5  RESANL 

Subroutine  RESANL  restarts  the  orbit  generator  through  a  call  to  subrou¬ 
tine  RESRNK,  a  Runge-Kutta  integrator,  and  through  initializing  interpolators. 
The  needed  incoming  parameters  are  input  through  RESANL’s  argument  listing 
or  COMMON  blocks.  They  are  described  in  Table  2.7.  All  of  RESANL’s  out¬ 
put  is  through  COMMON  blocks.  These  variables  arc  listed  in  Table  2.8  [55,15]. 
Subroutine  RESANL  is  called  by  BEG ANL  at  the  beginning  of  the  program  to  ini¬ 
tialize  the  integrator.  Subroutine  RESANL  can  also  be  called  by  the  main  program 
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Table  2.5.  ORBANL  Input 


Variable 

Definition 

OBSTIM 

Request  time  from  epoch,  A.l  sec 

/ANLPRM/ 

NSOLVE 

Number  of  solve  for  parameters 

/AVRHST/ 

DIRECT 

Direction  of  integration 
=  1,  forward 
=  -1,  backward 

/SPINTG/ 

INTPOS 

Interpolate  for  position/ velocity? 

=  1,  yes 
=  2,  no 

/SPREAL/ 

PVEND 

End  of  short-arc  interpolator  interval,  sec 

Table  2.6.  ORBANL  Output 


Variable 

Definition 

POS 

Position  vector  of  satellite 

VEL 

Velocity  vector  of  satellite 

OSCELM 

Osculating  equinoctial  elements  of  satellite 

AVRELM 

Mean  equinoctial  elements  of  satellite 

PVDRV 

Partial  derivatives  of  current  position/ velocity 
with  respect  to  solve  for  parameters 

AVRDRV 

Partial  derivatives  of  current  mean  equinoctial 
elements  with  respect  to  solve  for  parameters 

Does  external  event  stop  propagation? 

OBSTIM 

Request  time  from  epoch,  A.l  sec 

/SPREAL/ 

PVEND 

End  of  short-arc  interpolation  interval 
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Table  2.7.  RESANL  Input 


Variable 

Definition 

ELEMNT 

Commanded  new  mean 
equinoctial  orbital  elements 

PARDRV 

Partial  derivatives  of  the  mean  equinoctial 
elements  with  respect  to  solve  for  parameters 

RESTIM 

Time  elapsed  from  epoch 

/  AVRHST  / 

DIRECT 

Direction  of  integration 
=  1,  forward 
=  -1,  backward 

/MACHINE/ 

DBLNUL 

Null  value 

DBLMAX 

Maximum  value 

DBLMIN 

Minimum  positive  value 

Table  2.8.  RESANL  Output 


Variable 

Definition 

/SPREAL/ 

SPBEG 

Beginning  of  auxiliary  interpolator  interval 

SPEND 

End  of  auxiliary  interpolator  interval 

Beginning  of  position/velocity 
interpolator  interval 

PVEND 

End  of  position/velocity 
interpolator  interval 
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after  an  impulsive  burn  to  restart  the  integrator.  If  the  burn  model  is  not  impulsive, 
then  the  main  program  must  start  with  a  call  to  BEGANL  which  will  then  call 
RESANL  to  restart  the  integrator. 
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Chapter  3 

The  MEANELT  Program 


The  MEANELT  computer  program  is  an  existing  stationkeeping  tool  received 
from  Aerospace  Corporation.  It  provides  a  long  term  simulation  of  a  satellite’s  station¬ 
keeping  capability.  This  program  utilizes  the  Draper  Semianalytic  Satellite  Theory 
Standalone  Orbit  Propagator  (DSST  Propagator)  to  propagate  the  satellite  orbit  in 
time  until  the  satellite  needs  to  expend  an  impulsive  maintenance  burn.  Then,  the 
propagator  is  reinitialized  and  propagation  continues  until  the  stop  time  is  reached 
or  the  fuel  budget  is  depleted.  The  original  MEANELT  program  was  written  by  R. 
G.  Hopkins  of  the  Aerospace  Corporation  Astrodynamics  Department.  It  is  writ¬ 
ten  in  FORTRAN  77  in  51  modules  (3,  661  lines  of  code)  in  addition  to  the  DSST 
Standalone  subroutines  described  in  Chapter  2. 

3.1  Current  Capabilities 

The  stationkeeping  capabilities  of  MEANELT  include  maintaining  the  longitude  of 
perigee,  apogee,  or  ascending  node,  or  maintaining  the  semi-major  axis.  The  program 
ean  maintain  any  one  of  the  longitude  variables  within  a  time  period  of  up  to  ten 
years.  The  nominal  value  of  the  maintained  clement  is  allowed  to  change  over  time 
and  the  actual  maintained  longitude  will  be  a  linear  change  from  one  nominal  value 
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to  the  other.  This  program  is  robust  in  that  it  can  tolerate  stable  and  unstable  points 
of  an  orbit,  as  for  a  geosynchronous  satellite,  where  perturbing  accelerations  change 
sign.  There  are  no  assumptions  in  MEANELT  in  calculating  the  nodal  period  or  rate 
of  the  longitude  drift  [44]. 

MEANELT  also  has  the  ability  of  long  term  semi-passive  eccentricity  control.  This 
is  attained  through  varying  where  impulsive  changes  in  velocity  are  expended  and  the 
magnitude  of  the  impulsive  maneuver  to  maintain  the  longitude  of  perigee,  apogee, 
or  ascending  node,  or  the  semi-major  axis.  The  default  is  to  minimize  the  change 
in  velocity  for  longitude  maintenance,  but  there  is  also  a  capability  to  minimize  the 
change  in  eccentricity  or  to  keep  the  eccentricity  from  exceeding  a  nominal  value  of 
eccentricity,  eo  [45], 

3.2  Description  of  Original  Software 

Only  the  program  flow  for  maintaining  the  longitude  of  the  ascending  node  will 
be  discussed  since  that  capability  applies  to  maintenance  of  an  Earth  Observation 
satellite  in  a  sun  synchronous,  repeat  ground  track,  and  frozen  orbit.  By  maintaining 
the  longitude  of  the  ascending  node  of  every  M  satellite  orbits,  or  every  N  nodal 
days,  a  repeat  ground  track  orbit  can  be  maintained  using  the  program  MEANELT. 
The  total  complexity  of  MEANELT  is  much  too  great  for  the  scope  of  this  thesis. 

3.2.1  MAIN 

The  main  program  for  MEANELT  first  calls  two  subroutines,  SETDAF  and 
SETINP,  to  initialize  the  parameters  needed  by  the  orbit  propagator,  ie.  the  DSST 
Standalone,  and  the  stationkeeping  portion  of  the  program.  These  two  subroutines  are 
analogous  to  the  DSST  Standalone  subroutine  SETELM  described  in  Section  2.2.1. 
The  main  program  then  sets  variables  to  their  default  values  and  then  reads  nine 
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namelist  input  files  to  adjust  the  input  variables  of  the  default  case  to  the  desired 
case  simulation.  These  namelist  files  parallel  the  function  of  subroutine  INTANL 
of  the  DSST  Standalone.  The  nine  input  namelist  files  are  FRCSET,  AVRSET, 
SPSET,  DVRSET,  STKPSET,  MEANSET,  PROPSET,  OUTSET,  and  FU- 
ELSET.  The  namelist  FRCSET  specifies  the  geopotential,  spacecraft  physical,  solar 
radiation  pressure,  and  atmospheric  drag  constants.  It  also  includes  the  force  model 
used  for  the  orbital  perturbations.  AVRSET  contains  parameters  for  the  averaged 
equations  of  motion  including  propagation  error,  resonant  terms,  and  J2  gravitational 
potential  effects.  SPSET  indicates  the  parameters  that  determine  if  osculating  ele¬ 
ments  will  be  calculated  from  evaluated  short-periodic  terms.  DVRSET  determines 
if  partial  derivatives  of  the  current  mean  state  with  respect  to  the  epoch  mean  state, 
in  equinoctial  elements,  will  be  calculated.  STKPSET  sets  the  stationkeeping  pa¬ 
rameters.  MEANSET  defines  the  initial  mean  elements,  coordinate  system,  and 
reference  frame  for  the  repeat  orbit  and  stationkeeping  portion  of  the  program,  while 
the  namelist  PROPSET  defines  the  epoch  mean  elements,  coordinate  system,  ref¬ 
erence  frame,  and  output  format  for  the  propagation  portion  of  the  program,  ie.  the 
DSST  Standalone.  OUTSET  specifies  the  printed  and  plotted  output  to  be  pro¬ 
duced.  Lastly,  the  namelist  FUELSET  determines  if  a  fuel  budget  is  to  be  kept  and 
what  the  fuel  initial  conditions  are. 

The  main  MEANELT  program  finally  calls  subroutine  SYNKRO.  It  is  subroutine 
SYNKRO  that  acts  as  the  simulation  driver. 

3.2.2  SYNKRO 

Subroutine  SYNKRO  is  the  main  driver  for  the  orbit  maintenance  simulation. 
SYNKRO  first  converts  the  inputs  from  the  namelists  to  units  needed  by  MEANELT 
and  the  DSST  Standalone  Propagator.  This  driver  then  calls  two  separate  compo- 
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nents,  REPEAT  and  STKEEP  (see  Figure  3.1). 


Figure  3.1.  Top  Level  Driver 


Subroutine  REPEAT  adjusts  the  epoch  time  and  epoch  elements,  input  from  the 
namelists,  to  satisfy  user-imposed  repeatability  requirements  and  initial  perigee  alti¬ 
tude.  After  REPEAT  is  finished,  the  driver  SYNKRO  calls  subroutine  STKEEP 
to  coordinate  the  orbit  propagation,  print  the  output,  and  perform  stationkeeping  ma¬ 
neuvers  when  required.  Subroutine  STKEEP  is  the  real  ‘meat’  of  the  MEANELT 
stationkeeping  program  [44]. 

3.2.3  REPEAT 

Subroutine  REPEAT  adjusts  the  epoch  time,  epoch  elements,  and  input  from 
the  namelists  to  satisfy  user-imposed  repeatability  requirements  and  initial  perigee 
altitude.  To  determine  these  changes,  REPEAT  calls  subroutines  TNJECT  and 
WIEDER  (see  Figure  3.2).  Subroutine  TNJECT  initializes  the  DSST  Standalone 
and  adjusts  the  initial  epoch  time  to  achieve  the  desired  Earth  fixed  longitude  at  the 
ascending  node.  It  is  subroutine  TNJECT  that  calls  the  DSST  Standalone  interface 
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TNJECT 

WIEDER  | 

Figure  3.2.  Subroutine  REPEAT 


subroutine  BEGANL.  A  subroutine  with  a  shadowed  box  in  the  flow  diagram  signi¬ 
fies  a  probable  site  of  major  modifications  in  MEANELT  for  a  simulation  similar  to 
LANDSAT  6.  If  the  desired  orbit  is  not  only  constrained  by  a  repeat  ground  track 
but  also  by  sun  synchronous  and  frozen  orbit  constraints,  then  adjustment  of  the 
initial  epoch  time  in  TNJECT  is  not  needed  since  the  epoch  time  given  will  already 
achieve  the  needed  orbital  elements.  Subroutine  WIEDER  should  be  stubbed  out 
all  together  for  the  same  reasons.  Using  WIEDER  iteratively,  the  semi- major  axis 
can  be  adjusted  to  fit  the  repeatability  requirements  with  a  Newton- Raphson  tech¬ 
nique.  If  an  initial  value  of  perigee  altitude  is  specified  by  the  user,  further  iteration  in 
WIEDER  can  be  performed  to  adjust  both  the  semi-major  axis  and  eccentricity  [44], 

3.2.4  STKEEP 

Subroutine  STKEEP  coordinates  the  orbit  propagation,  prints  the  output, 
and  performs  maintenance  maneuvers  when  required.  This  is  the  real  ‘meat’  of  the 
MEANELT  stationkeeping  program.  To  accomplish  these  functions,  STKEEP  calls 
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subroutines  APOELM,  OUTPR,  BDLAM,  and  COREKT,  refer  to  Figure  3.3. 
Subroutine  APOELM  uses  Lagrangian  interpolation  in  the  z-coordinate  to  deter- 


Figure  3.3.  Subroutine  STKEEP 


mine  the  time  of  the  ascending  node  crossing.  Inside  of  this  subroutine,  the  DSST 
Standalone  determines  the  mean  orbital  elements  at  that  node  crossing  time  through 
a  call  to  the  DSST  Standalone  ORBANL.  The  subroutine  APOELM  also  makes 
a  call  to  the  DSST  Standalone  subroutine  RESANL  if  an  impulsive  maneuver  was 
just  completed.  Subroutine  OUTPR  outputs  the  data  concerning  orbit  propagation 
and  stationkeeping  maneuvers  to  a  file  and  writes  various  data  to  arrays  for  later  plot¬ 
ting.  The  subroutine  that  compares  the  geographic  longitude  of  the  ascending  node 
to  the  range  of  allowed  values  is  BDLAM.  If  no  stationkeeping  is  needed,  propaga¬ 
tion  continues.  Subroutine  COREKT  performs  the  actual  stationkeeping  maneuver 
if  needed  and  maintains  the  fuel  consumption  of  the  spacecraft.  COREKT  along 
with  WIEDER  iteratively  determine  the  value  of  the  semi-major  axis  and  eccen¬ 
tricity  to  produce  the  desired  variation  in  the  longitude  of  the  ascending  node.  Here 
again,  subroutine  WIEDER  could  be  stubbed  out  and  a  more  direct  calculation  of 


74 


determining  the  new  semi- major  axis  substituted  [44]. 

It  is  in  subroutine  COREKT  that  the  long  term  passive  eccentricity  control  takes 
place  if  desired.  COREKT  computes  the  magnitude  of  the  maintenance  burn  needed 
to  adjust  the  semi-major  axis  to  maintain  the  repeat  ground  track,  ie.  the  longitude 
of  the  ascending  node.  This  burn  may  occur  at  many  places  in  the  orbit  to  achieve 
the  change  in  the  semi-major  axis  desired.  Through  varying  where  this  impulsive 
maneuver  takes  place  and  varying  the  magnitude  of  the  burn,  subroutine  COREKT 
is  able  to  either  minimize  the  velocity  expended  for  longitude  maintenance,  minimize 
the  change  in  eccentricity,  or  keep  the  eccentricity  from  exceeding  a  nominal  value  of 
e0.  In  the  default  case,  IALT  =  0,  where  the  change  in  velocity  is  minimized,  the 
change  in  eccentricity  is  approximately: 

Ae  «  (1  —  e)A a/a  (3.1) 

In  this  case,  the  change  in  velocity  is  expended  at  perigee  to  adjust  the  orbit  altitude 
and  also  affect  the  eccentricity.  This  type  of  maneuver,  that  changes  both  the  semi¬ 
major  axis  and  the  eccentricity,  is  not  suitable  for  an  orbit  with  a  small  eccentricity 
when  the  change  in  the  semi-major  axis  due  to  perturbations  is  less  than  zero,  ie.  drag 
is  the  dominant  perturbation.  In  this  case  where  the  semi-major  axis  is  decreasing,  the 
burn  tries  to  force  the  eccentricity  to  a  value  less  than  zero.  Therefore,  the  following 
two  capabilities  were  added  to  accommodate  orbits  with  small  eccentricities.  If  the 
change  in  eccentricity  during  a  burn  is  to  be  minimized,  IALT  =  1,  the  semi-major 
axis  will  be  changed  approximately  independently  of  the  other  orbital  elements,  see 
Table  1.6.  Half  of  the  AV  for  the  altitude  adjust  will  be  spent  at  perigee  and  half  at 
apogee  for  small  A  a/ a.  The  change  in  eccentricity  for  this  case  is  approximately: 

Ae  %  — e(A a/a)  (3.2) 

This  alternate  burn  will  cause  only  slightly  more  fuel  to  be  used  than  in  the  default 
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option.  If  the  eccentricity  is  to  be  kept  from  exceeding  a  set  nominal  value  of  eccen¬ 
tricity,  IALT  =  2,  a  combination  of  the  above  two  cases  will  be  used.  If  the  current 
eccentricity  is  less  than  the  nominal  value,  then  the  minimum  change  in  eccentricity 
burn  will  be  simulated.  If  the  eccentricity  is  greater  than  the  nominal  value,  then  the 
change  in  velocity  will  be  expended  at  perigee  if  Aa  <  0  or  at  apogee  if  A  a  >  0.  The 
latter  will  lead  to  a  change  in  eccentricity  of: 

Ae  ~  —  (1  +  e)Aa/a  (3.3) 

This  third  option  will  succeed  if  the  needed  change  in  the  semi-major  axis  is  small 
or  maneuvers  are  made  infrequently.  In  addition,  maintaining  the  eccentricity  below 
a  nominal  value  leads  to  a  substantial  increase  in  the  velocity  expended  and  the 
simulation  computational  time  needed  if  the  change  in  the  semi-major  axis  due  to 
orbit  perturbations  is  less  than  zero. 

Subroutine  COREKT  assumes  an  impulsive  burn  model  which  may  be  replaced 
by  actual  burn  models  dependent  on  time  that  are  discussed  further  in  Section  4.2.  In 
addition,  the  original  repeatability  requirements  can  be  replaced  with  Bruce  Baxter’s 
algorithms,  see  Section  4.1.2  for  a  specific  mathematical  discussion  of  his  concepts. 

After  the  new  orbital  elements  and  required  impulsive  Av  are  calculated,  orbit 
propagation  continues.  This  is  accomplished  by  STKEEP  making  the  necessary  call 
to  APOELM  to  restart  the  orbit  propagator  through  subroutine  RESANL. 
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Chapter  4 

Maneuver  Planning  Algorithms 


4.1  Orbit  Design  Concepts 


Earth  Observation  satellites  use  many  common  orbit  design  techniques  to  improve 
their  instrumentation  capabilities.  Three  of  the  most  common  orbit  design  concepts, 
sun  synchronous  orbit,  repeat  ground  track,  and  frozen  orbit,  were  introduced  in 
Chapter  1  and  are  discussed  in  more  mathematical  detail  here. 

4.1.1  Sun  Synchronous  Orbit 


A  sun  synchronous  orbit  keeps  the  orbit  plane  at  a  constant  angle  to  the  sun  by 
using  the  geopotential  perturbations  of  the  Earth  to  force  the  line  of  nodes  westward 
along  the  equator  approximately  one  degree  per  day.  The  line  of  nodes  is  defined  as 
the  intersection  of  the  orbit  plane  and  the  Earth’s  equator.  This  line  of  nodes  is  often 
stated  as  the  nearly  constant  local  time  of  ascending  node  crossing  rather  than  an 
angular  measurement  of  longitude.  The  o’clock  angle  is  the  nearly  constant  angular 
value  between  the  line  of  nodes  and  the  Earth/sun  line. 

Through  the  application  of  the  Earth’s  geopotential  to  Lagrange’s  Planetary  Equa¬ 
tions,  the  mean  nodal  rate  is  calculated  [18]: 


dtt 

dt 


RP 


,«(!  ~e2) 


(4.1) 
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where  Cl  is  the  mean  longitude  of  the  ascending  node,  the  angle  of  the  line  of  nodes 
from  the  vernal  equinox.  The  variable  J2  is  the  harmonic  coefficient  of  the  Earth’s 
oblateness  equal  to  1.08263xl0-3  and  the  gravitational  parameter  of  the  Earth,  /z,  is 
equal  to  398,601.2  km3/  sec2.  The  mean  equatorial  radius  of  the  Earth,  denoted  Re, 
is  equal  to  6378.135  km.  These  values  are  approximate  and  they  depend  on  the  actual 
gravity  model  used.  The  semi-major  axis,  a,  the  inclination,  z,  and  the  eccentricity, 
e,  of  the  orbit  are  variables  in  the  common  Keplerian  orbital  element  set.  By  setting 
the  mean  nodal  rate  to  sun  synchronous: 

dCl  360°  7 r  year  day  hr 

dt  year  180°  365.25 day  24 hr  3600sec 

=  1.991021  x  10-7rad/sec  (4.2) 

and  assuming  that  the  orbit  is  near  circular,  ie.  e  ss  0,  Equation  4.1  may  be  altered  to 
solve  for  the  sun  synchronous  relationship  between  the  inclination  and  the  semi-major 

axis: 


i  =  arccos  (-4.773621  x  10_1V/2)  (4.3) 

The  graph  of  this  relationship  is  shown  in  Figure  4.1. 

One  of  the  primary  perturbations  not  considered  in  the  design  of  the  sun  syn¬ 
chronous  orbit  is  the  third  body  effects  due  tc  the  point  mass  of  the  sun.  The  effect 
of  solar  gravity  on  a  sun  synchronous  spacecraft  causes  the  orbital  inclination  to 
increase  or  decrease  depending  on  whether  the  ascending  node  is  PM  or  AM,  respec¬ 
tively.  This  change  in  the  inclination  creates  a  change  in  the  nodal  precession  rate 
(Equation  4.1)  which  causes  a  shift  in  the  expected  ground  track.  This  can  cause  an 
increase  in  ground  track  maintenance  if  the  ground  track  is  to  repeat  [26]. 
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Earth's  Surface 


Figure  4.1.  Sun  Synchronous  Orbit:  a  vs.  i 

4.1.2  Repeat  Ground  Track 

A  satellite  ground  track  is  defined  as  the  locus  of  points  traced  out  on  the 
Earth’s  surface  directly  beneath  the  spacecraft  orbit.  This  ground  track  does  not 
noimally  repeat  itself,  ie.  retrace  a  previous  orbit’s  ground  track,  by  reason  of  the 
Earth  rotating  underneath  the  satellite  orbit  while  the  line  of  nodes  rotates  around 
the  Earth’s  equator  due  to  perturbations.  A  commensurability  condition  for  a  repeat 
ground  track  may  be  determined  by  using  an  equation  from  Baxter  [5]: 

2-kN  =  M(ue~  ti0)PNo  (4.4) 

where  N  is  the  number  of  nodal  days  until  the  ground  track  repeats,  M  is  the  number 
of  orbits  in  the  repeat  cycle,  uie  is  the  approximately  constant  rotation  rate  of  the 
Earth,  f20  is  the  initial  mean  rotation  rate  of  the  line  of  nodes,  and  P^0  is  the  initial 
mean  nodal  period  of  the  satellite.  A  nodal  day  is  defined  as  the  time  for  the  Earth 
to  rotate  360°  with  respect  to  the  line  of  nodes.  A  solar  day  is  defined  as  the  time 
for  the  Earth  to  rotate  360°  with  respect  to  the  Earth-Sun  line.  If  the  satellite  orbit 
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is  sun  synchronous,  a  solar  day  is  approximately  the  same  as  a  nodal  day  since  the 
line  of  nodes  is  at  a  constant  angle  to  the  sun.  If  the  satellite  is  near  circular,  ie. 
eccentricity  is  approximately  zero,  the  equation  for  the  nodal  period  is  [78,18]: 

PN  =  Pk  1  ~ (~)  (6  —  5 sin2  z)  +  0(e2)  (4.5) 

where  Pk  is  defined  as  the  Keplerian  (two  body)  orbital  period  [4]: 

Pk  =  27'^[j  (4'6) 

where  a  is  the  semi-major  axis  and  i  is  the  inclination  and  both  values  are  taken  at  the 
node  crossing.  As  previously  stated  in  Equation  4.1,  the  variable  J2  is  the  harmonic 
coefficient  of  the  Earth’s  oblateness  equal  to  1.08263xl0-3,  the  mean  equatorial  radius 
of  the  Earth,  denoted  Re ,  is  equal  to  6378.135  km,  and  the  Earth’s  gravitational 
parameter,  /i,  is  398601.2  km3/sec2.  Again,  these  values  are  approximate  and  they 
depend  on  the  actual  gravity  model  used. 

By  observing  the  commensurability  condition  for  a  repeating  ground  track  in 
Equation  4.4,  the  variable  u>e  represents  the  eastward  drift  of  the  Earth  while  the 
variable  represents  the  natural  westward  drift  of  the  line  of  nodes  from  orbital 
perturbations,  primarily  J2.  The  cumulative  drift  difference  must  be  some  integer 
multiple  of  2n  after  M  revolutions  for  the  ground  track  to  exactly  repeat  itself.  By 
manipulating  Equation  4.4,  an  equation  for  the  longitudinal  shift  of  the  line  of  nodes 
per  orbitn'  revolution,  5(0),  is  calculated  from  the  mean  nodal  period  and  nodal  rate 
at  epoch  [5]. 

.  N  27 r 

5(0)  =  (we  -  )Pn0  =  TvT  -  g-  (4-7) 

where  Q  is  the  number  of  revolutions  per  day  of  the  satellite.  An  illustration  of 
this  longitudinal  node  drift  with  a  non-rotating  Earth  can  be  viewed  in  Figure  4.2. 
The  variable  N  in  this  coordinate  system  can  now  be  viewed  as  the  number  of  times 
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Reu  0 

Figure  4.2.  Node  Drift 
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that  the  line  of  nodes  will  drift  around  the  equator  of  a  stationary  Earth  until  it 
returns  to  its  original  position.  If  the  variable  S(j)  is  the  true  longitude  drift  per 
orbital  revolution  for  the  jth  revolution,  then  the  value  DS(j )  is  the  instantaneous 
differential  drift  in  longitude  from  the  original  longitudinal  drift  defined  by: 

DS(j)  =  S(j)-S(  0) 

—  "b )  PN] 

=  ue(PNj  ~  P N0)  —  foj Pn}  +  Pn0  (4.8) 

By  assuming  that  the  longitudinal  rate  is  linear  over  the  time  of  one  orbit  revolu¬ 
tion,  the  nodal  rate  can  be  equated: 

n,  =  ,4.9) 

Now,  substitute  Equation  4.9  into  Equation  4.8: 

DSj  =  ue(PNj  -  PNo)  -  [(fy  -  a,_i)  -  n0PNo)  (4.10) 

From  this  equation,  the  value  DSj  can  be  easily  determined  from  the  satellite  ephemeris 
file.  An  ephemeris  is  a  log  of  the  past  positions  of  the  satellite  in  orbit.  Given  the 
longitude  of  the  ascending  node,  the  time  of  crossing  can  be  interpolated  from  the 
ephemeris  file.  Using  this  time,  the  whole  satellite  state  can  be  found  from  the 
ephemeris  file.  From  this  satellite  state,  the  nodal  period  and  nodal  rate  can  be 
calculated,  thus  the  value  DSj  is  determined.  By  letting: 

SPn  =  I\-Pn0  (4.11) 

5Q  =  (Qj  -Uj-i)  -SWo  (4.12) 

the  instantaneous  differential  drift  in  longitude  of  the  ascending  node,  DS,  is  now 
equated: 

DS  =  uebPN  -  bil  (4.13) 
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The  variable  DS  can  also  be  viewed  as  a  derivative,  therefore: 


S  =  So  +  lN (DS)dR 
Jo 


(4.14) 


where  S  is  the  cumulative  drift,  R  is  the  revolution  number  and  Sq  is  the  initial 
displacement  of  the  longitude  of  the  ascending  node  with  respect  to  the  nominal 
value.  This  value,  Sq,  is  calculated: 


N 

So  =  J2S(0)  =  NS(0)  (4.15) 

»=i 

The  variable  DS  is  approximately  a  linear  function  of  time,  ie.  the  rev  number, 
because  it  is  primarily  a  function  of  the  value  6Pn  which  is  an  almost  linear  function 
of  the  slowly  changing  semi-major  axis  principally  due  to  atmospheric  drag.  Since 
the  change  in  the  longitude  of  the  ascending  node  is  small  compared  to  the  value  for 
the  longitude  of  the  ascending  node,  the  value  DS  can  expressed  linearly: 


DS  —  Qq  +  ol\R 


(4.16) 


where  the  a  coefficients  are  real  numbers.  When  this  linear  model  is  substituted  into 
Equation  4.14  and  integrated,  the  cumulative  drift  is  now  a  quadratic: 

S  =  S0  T  a\R  +  a2R2 

—  Qq  T  0,1  R  +  (I2R?  (4.17) 


The  values  for  cumulative  drift,  S ,  may  be  smoothed  with  a  sequential  Least  Squares 
filter  to  remove  the  medium  and  short  periodics  and  to  calculate  the  coefficients: 
ao, oi,a2-  Now,  the  equation  for  DS  is  solved  with  q0  =  ai  and  a i  =  2a2- 

If  Srcf  is  the  reference  longitude  about  which  the  ground  track  is  to  be  maintained, 
then  the  values  Sref±S*  are  the  upper  and  lower  boundary  longitudes  where  S*  is  the 
ground  track  boundary  region  around  Sref  ■  Equation  4.17  gives  the  station  keeping 
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ability  within  Sre/  ±5*,  or  Sup  andSiow  Since  atmospheric  drag  decreases  the  semi¬ 
major  axis,  a  ground  track  correction  will  always  be  on  the  lower  boundary.  For 
example,  by  allowing  S  —  S[ow,  Equation  4.17  may  be  solved  for  the  number  of  orbit 
revolutions  for  the  ground  track  to  drift  from  the  lower  boundary,  reach  the  maximum 
drift,  and  return  to  the  lower  boundary,  ie.  let  R  =  Riow 


Slow  flo  d"  QlRlow  +  a2Rlow 


d2 Rlow  QlRlow  d~  do  Slow  0 


di  i  yd'f  4a2(d0  Sjow) 

ftiow  —  z  loj 

^0*2 

To  find  the  number  of  orbit  revolutions  until  the  maximum  or  minimum  of  the 
cumulative  drift,  S,  the  partial  derivative  of  S  with  respect  to  R  is  set  equal  to  zero. 

ss 

6R 

ai  1 2d2 R  ~  0 

-a. 


R  = 


202 


(4.19) 


Now  the  cumulative  drift  extremum,  5e,  is  determined  by  substituting  Equation  4.19 
into  Equation  4.17: 


c  ,  ai(-oi)  .  02 a? 

s‘  ~  a°  +  ^r“  +  i3' 


—  do  — 


4do 


(4.20) 


Since  the  instantaneous  drift  rate  is  the  slope  of  S  at  that  point,  the  drift  rate  at 
Riow  where  the  orbit  correction  is  to  take  place  will  equal: 


6S  , 

—  =  a{  +  2  a2Riow 

!  Ii—RUnu 


(4.21) 


This  drift  rate  must  not  only  be  stopped  but  reversed  so  that  the  drift  will  approach 
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the  extrema  value,  Se ,  again.  Therefore,  the  drift  correction,  DC,  can  be  calculated: 


DC  =  “2  1  H 1  =  ~2(ai  +  2 a2Rlow)  (4.22) 

\6RJr^r1ow 

Now,  the  drift  correction  must  be  related  to  the  orbital  elements  so  that  the  amount  of 
Ah'  may  be  calculated  to  maintain  the  ground  track.  This  relationship  is  furnished  by 
Equation  4.8  since  it  defines  the  instantaneous  drift  rate  in  terms  of  the  perturbations 
in  the  nodal  period,  6Pn,  and  the  value  6CI.  A  first  order  Least  Squares  fit  to: 

6D  =  b0PblR  (4.23) 


using  ephemeris  data  with  Equation  4.12  to  attain  values  for  6Q  will  remove  the  linear 
trend  caused  by  6Q  from  the  quadratic  drift  correction  equation.  This  is  possible  since 
a  small  change  in  the  nodal  period  will  have  a  negligible  effect  of  6 Q.  To  induce  the 
necessary  drift  rate  set  DS  =  DC: 


DC  =  DS  =  ojJPn  -  sn 
DC  =  uje5P n  —  b0  —  b\Riow 

6Pn 


— 2(<2j  T  2a2Riow) 

— 2(oj  +  2a2Riou)) 

bo  —  2oi  —  (4ci2  —  bi)Riow 


(4.24) 


From  Equation  4.5,  the  nodal  period  is  primarily  a  function  of  the  semi-major  axis, 
and  to  a  lesser  order  the  inclination  and  eccentricity.  Since  the  nodal  period  is  basi¬ 
cally  a  function  of  the  semi-major  axis,  especially  for  a  sun  synchronous  orbit  where 
the  inclination  is  near  polar,  a  valid  assumption  to  be  made  is  that  the  nodal  period 
can  be  approximated  by  the  Kcplerian  period,  ie.  Ps  ~  Pk  where  Pk  is  given  in 
Equation  4.6.  Differentiating  the  Keplerian  period  determines  6Pu- 
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(4.26) 


Recalling  the  vis-viva  integral: 


v 


2 


(4.27) 


By  assuming  a  circular  orbit  where  the  radial  distance,  r,  equals  the  semi-major  axis, 
a,  the  vis  viva  integral  is  now: 


v2  =  £ 


(4.28) 


Now  differentiate  and  rearrange  Equation  4.28: 

2v8v  =  —  ■J'—Sa 

8a  a  .  8v 

—  =  —2 —v8v  =  2 — 
a  n  v 


(4.29) 


Now  substitute  this  definition  for  ~  into  Equation  4.25  and  solve  for  AV: 


8Pn 

6v 

AV  =  |<H 


\Pn  f — 2— 

2  V  v 


v  8  Pm 

'3^ 
v 


'■iPpjtUf 


[6q  —  2ai  —  (4a2  —  &i)/?iOU;] 


(4.30) 


This  is  an  approximate  value  for  the  magnitude  of  AV  needed  to  correct  the  semi¬ 
major  axis  since  the  nodal  period  was  approximated  with  the  Keplerian  period. 

From  Equation  4.17,  the  parabolic  relationship  between  the  orbit  revolutions  and 
the  amount  of  longitudinal  drift  can  be  graphed,  thus  showing  the  orbit  correction 

cycle. 

From  this  illustration,  it  is  easy  to  see  that  the  decrease  in  the  semi-major  axis 
from  atmospheric  drag  causes  the  ground  track  adjust  to  always  occur  at  the  lower 
boundary.  As  the  AV  produces  a  positive  8 Pm,  the  line  of  nodes  drifts  west,  but  as 
PN  decreases  due  to  drag  the  orbit  will  reach  a  point  where  PN  =  Pc  at  the  extremum 
and  the  line  of  nodes  will  then  drift  east  until  the  orbit  is  corrected  again. 
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Ground  Track  Drift  Longitude 


Orbit  Revolutions 


Figure  4.3.  Orbit  Correction  Cycle 
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It  would  be  an  optimal  burn  if  the  maximal  cumulative  drift,  Se,  was  the  same 
as  the  upper  boundary,  Sup.  But  because  of  the  computational  lag  for  the  data  fit 
in  the  beginning  of  the  parabolic  arc  and  uncertainties  in  the  future  atmospheric 
drag,  the  burn  needs  to  under-correct  a  bit  to  avoid  having  to  correct  again  at  the 
upper  bound,  usually  needing  a  180°  yaw  or  pitch  attitude  maneuver.  Therefore,  an 
absolute  calculation  of  AV  is  impossible  [5]. 

There  are  two  types  of  errors  that  can  destroy  the  repeating  ground  track.  Pri¬ 
marily,  there  is  atmospheric  drag  that  decreases  the  semi-major  axis  and  causes  the 
ground  track  drift  trend  to  be  eastward.  An  illustration  of  this  drift  can  be  viewed 
in  Figure  4.4  [18]. 


Reference  Ground  Track 


Figure  4.4.  Ground  Track  Drift  due  to  Semi-major  Axis  Decay 
The  second  type  of  ground  track  drift,  due  to  changes  in  inclination,  affects  the 
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nodal  period  and  ground  track  errors  at  higher  latitudes.  A  typical  AM  sun  syn¬ 
chronous  orbit  will  have  a  slow  decrease  in  inclination  which  will  cause  a  slight  de¬ 
crease  in  the  rate  of  node  regression  and  a  slight  increase  in  the  nodal  period  [48]. 
Inclination  drift  will  also  cause  ground  track  error  at  higher  latitudes.  The  orbit  incli¬ 
nation  is  perturbed  mainly  by  solar  third  body  gravitational  effects,  thus  slanting  the 
orbit  more  and  upsetting  the  repeat  ground  track  at  the  upper  latitudes.  An  example 
of  this  type  of  ground  track  drift  can  be  viewed  in  Figure  4.5  [18].  An  illustration 


Reference  Ground  Track 


Figure  4.5.  Ground  Track  Drift  due  to  Inclination  Decay 

of  the  cumulative  effect  of  the  two  types  of  ground  track  errors  can  be  viewed  in 
Figure  4.6. 

The  semi-major  axis  is  periodically  boosted  with  a  burn  in  the  posigrade  direction 
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Reference  Ground  Track 


Figure  4.6.  Cumulative  Ground  Track  Drift 
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using  the  AV  determined  previously  is  Equation  4.30.  The  inclination  is  boosted 
with  an  out-of-plane  maneuver,  but  this  type  of  boost  is  done  less  frequently  than 
the  semi-major  axis  adjust.  As  mentioned  in  Section  4.1.1,  the  inclination  drift  in  a 
sun  synchronous  orbit  not  only  needs  to  be  adjusted  to  maintain  the  sun  synchronous 
condition  but  it  also  needs  to  be  adjusted  to  maintain  the  ground  track  at  the  upper 
latitudes.  The  out-of-plane  maneuver  needed  for  this  inclination  adjust  may  be  costly 
because  of  the  attitude  maneuver  necessary  to  place  the  thrusters  in  the  correct 
position.  For  a  sun  synchronous  orbit,  the  satellite  may  have  its  inclination  biased 
such  that  the  inclination  will  drift  towards  and  then  away  from  the  desired  inclination 
value,  thus  elongating  the  time  to  an  inclination  adjust  maneuver.  The  inclination 
change  per  day  of  a  sun  synchronous  satellite  can  be  estimated  [26]: 

/  =  +  cos  is)2  sin i0  sin  2Qoc  (4-31) 

day  8 n2  rj 

where  Q  is  again  the  number  of  orbits  per  day,  n  is  the  satellite  mean  motion  equal 
to  is  the  solar  gravitational  parameter,  rs  is  the  mean  distance  to  the  center 

of  the  sun  from  the  Earth’s  center,  i0  is  the  satellite  inclination  at  epoch,  ttoc  is 
the  o’clock  angle,  and  is  is  the  solar  obliquity  angle  approximately  equal  to  23.45°. 
This  change  in  inclination  is  very  slow,  thus  it  is  reasonable  to  make  a  first  order 
approximation  of  the  inclination: 

*  =  *0  +  IT  (4.32) 

Thus  the  cosine  of  the  inclination  can  be  calculated: 

cos  i  =  cos  io  cos  IT  —  sin  io  sin  IT 

«  cos  i0  —  IT  sin  i0  (4.33) 

The  latter  equation  comes  from  the  assumption  that  the  time  multiplied  by  the  rate 
of  change  of  the  inclination  in  radians  is  very  small.  Now  substituting  Equation  4.33 


91 


into  the  equation  for  the  nodal  rate,  Equation  4. 1  and  the  assumption  that  the  orbit 
is  near  circular: 


1  COS  2 


=  —  -J2  (^~^j  n[cosi0  -  IT sini0] 

=  ^0  +  2^2  (“J  n  IT  sin  i0 

If  it  is  assumed  that  the  mean  node  rate  is  a  linear  approximation: 


(4.34) 


—  Qq  T  QT 


(4.35) 


then  the  acceleration  of  the  line  of  nodes  is  calculated: 


HMD 


nl  sini0 


(4.36) 


By  integrating  Equation  4.35,  an  equation  for  the  line  of  nodes  can  be  found: 

n  -  Qo  +  tioT  +  ^-T2 

n  =  nss  +  6n  =  n0  +  {hss  +  6Ciq)t  +  ^ j2  niT2  sin i0  (4.37) 

where  Clss  is  the  node  regression  rate  for  a  sun  synchronous  orbit,  <5f2o  is  the  initial 
error  in  the  node  regression  rate,  and  Qss  is  the  instantaneous  line  of  nodes  for  a  sun 
synchronous  orbit  which  equals  f2o  f  f)ssT.  The  initial  error  in  the  line  of  nodes, 
is  calculated: 

6fl  =  6h0T  +  ^T2  (4.38) 

where  the  first  term  in  this  equation  is  due  to  launch  vehicle  and  orbit  injection 
errors  and  the  second  term  is  due  to  solar  gravitational  effects.  The  mairfenance 
requirements  for  the  inclination  drift  can  be  reduced  by  injecting  the  satellite  into  an 
orbit  with  an  inclination  bias  of  an  opposite  polarity  of  inclination  change  than  that 
caused  by  solar  gravity  [26], 
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4.1.3  Frozen  Orbit  Concept 


The  frozen  orbit  concept  is  based  on  minimizing  the  mean  eccentricity  and  ar¬ 
gument  of  perigee  rates  of  change  due  to  the  Eaxth’s  geopotential  perturbations.  A 
frozen  orbit,  where  the  eccentricity  and  argument  of  perigee  remain  approximately 
stable,  establishes  a  balance  between  the  secular  and  long  periodic  contributions 
caused  by  the  even  zonal  harmonics  and  the  long  periodic  contributions  caused  by 
the  odd  zonal  harmonics.  The  primary  harmonics  that  affect  the  rates  of  change  for 
the  eccentricity,  e,  and  argument  of  perigee,  u,  are  the  J2  and  J3  zonal  harmonic 
terms  [61].  In  the  orbit  analysis  for  SEASAT,  the  mean  element  rates  for  the  eccen¬ 
tricity  and  argument  of  perigee  arc  calculated  [21]: 


de 

3nRlJ3  sin  i 

dt 

2a3(l  —  e2)2 

du 

3  nJ2Rl  / 

dt 

a2(l  -  e2)  V 

('-  jsin2‘) 
-  ^  sin2  i  j  0 


cosu; 


0  =  1  + 


J$Re 


'sin2  i  —  e  cos2  sin  a; 


2J2a{l  -  e2) 


sinr 


(4.39) 

(4.40) 

(4.41) 


where  n  is  the  mean  motion  of  the  satellite  and  Re  is  the  Earth’s  mean  equatorial 
radius.  The  elements  in  these  equations  are  assumed  to  be  mean  elements.  Inspecting 
Equation  4.39,  the  mean  eccentricity  rate  will  approach  zero  for  an  equatorial  orbit  (ie. 
i  =  0, 180°),  an  orbit  at  critical  inclination  (ie.  icrit  =  63.4°,  116.6°),  or  an  orbit  which 
has  an  argument  of  perigee,  u,  at  ±90°.  If  the  orbit  inclination  is  constrained,  as 
for  sun  synchronous  orbits  or  Earth  coverage  constraints,  the  only  variable  available 
to  minimize  the  mean  eccentricity  rate  is  the  argument  of  perigee.  By  inspecting 
Equation  4.40,  the  mean  argument  of  perigee  rate  vanishes  when  the  orbit  is  at 
critical  inclination  or  when  the  variable  0  goes  to  zero.  For  the  same  reasons  as  in 
the  eccentricity,  only  ©  is  available  for  orbit  selection.  By  setting  0  to  zero,  the 
nominal  eccentricity,  eo,  that  will  minimize  the  rate  of  argument  of  perigee  can  be 
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determined: 

it  R 

e0  sa  —  — - - -  sin i  sin cj  +  O(e)  (4.42) 

2j2  a 

If  the  orbit  is  not  near  the  critical  inclination,  this  equation  is  mainly  driven  by  the 
J-2  term,  even  when  higher  zonal  harmonics  are  added.  This  means  that  the  nominal 
eccentricity  will  be  nearly  equal  to  l^|  ~  J2  %  .001.  Although  the  eccentricity  does 
change  when  higher  zonals  are  added  to  the  gravity  model,  the  eccentricity  never 
changes  more  than  twenty  percent.  If  the  orbital  inclination  is  near  critical,  the 
eccentricity  is  no  longer  approximately  equal  to  J2  and  can  reach  values  much  larger. 

Other  perturbations,  such  as  drag,  cause  the  eccentricity  to  drifi  from  thus  nominal 
value.  Therefore,  the  eccentricity  must  be  periodically  adjusted  with  maintenance 
burns.  In  orbit  analysis  for  NROSS,  it  was  determined  that  the  change  in  eccentricity 
would  be  very  small  over  time  from  perturbations  other  than  the  Earth’s  geopotential, 
and  this  small  change  could  be  easily  counteracted  during  the  required  semi-major 
axis  adjustment  burns.  An  additional  advantage  to  the  frozen  orbit  concept,  besides 
passive  control  of  the  eccentricity  and  argument  of  perigee,  is  that  all  global  and  local 
ai  ;  i I, tide  variations  are  minimized.  The  global  altitude  variation  is  minimized  through 
the  damping  of  the  long  period  motion  of  the  eccentricity.  The  local  altitude  variation 
is  minimized  since  the  argument  of  perigee  oscillates  around  its  stable  value  of  ±90° 
instead  of  rotating  about  the  orbit  [61]. 

4.2  Maneuver  Models 

4.2.1  Impulsive  Targeting  Model 

Ground  track  drift  necessitates  thruster  burns  to  maintain  the  perturbed  orbital 
elements,  specifically  the  semi-major  axis,  o,  and  the  orbital  inclination,  i.  The 
thruster  models  used  in  the  LANDS/vT  6  mission  profile  are  not  modelled  as  impulsive 


burns,  but  an  impulsive  model  can  be  used  initialize  the  burn  parameters  of  the  burn 
model.  These  parameters  can  be  refined  in  numerical  or  analytical  integration  of  the 
equations  of  motion  for  a  continuous  thrust  model. 

Inclination  Maintenance 

As  discussed  above,  the  inclination  is  mainly  perturbed  in  a  sun  synchronous  orbit 
by  solar  gravity.  This  change  in  orbit  inclination  not  only  causes  displacement  of  the 
ground  track  at  higher  latitudes,  but  it  also  affects  the  nodal  period,  see  Equation  4.5. 
The  inclination  burns  to  correct  the  perturbations  are  expected  to  take  place  every 
6  months  if  the  ground  track  is  maintained  within  ±1  km,  or  every  30  months  if  the 
ground  track  is  maintained  within  ±5  km  [62].  Since  the  LANDSAT  6  spacecraft 
will  have  thrusters  in  only  one  direction,  an  inclination  adjust  maneuver  will  require 
an  approximately  90°  yaw  to  position  the  thrusters  in  the  correct  direction.  The 
change  in  inclination,  i,  to  first  order  in  the  eccentricity  can  be  determined  from  the 
equation  [62]: 

A i  =  —Avn  cos  u  +  0(e2)  (4-43) 

na 

where  n  is  the  mean  motion  of  the  satellite,  Avn  is  the  magnitude  of  the  change  in 
velocity  in  the  positive  normal  direction,  and  u  is  the  argument  of  latitude  which 
equals  the  sum  of  the  true  anomaly,  /,  and  the  argument  of  perigee,  u.  This  equation 
can  be  inverted  to  determine  the  impulsive  value  of  the  velocity  needed  in  the  normal 
direction  to  achieve  the  desired  change  in  the  inclination  [62]: 

Tld 

Avn  = - Ai  (4.44) 

cos  u 

According  to  McClain  [62],  this  yaw  angle  to  position  the  thrusters  for  an  incli¬ 
nation  burn  is  not  exactly  90°.  Since  an  exact  out  of  plane  burn  would  result  in 
a  change  in  the  semi-major  axis,  an  unnecessary  drift  in  the  ground  track  could  be 
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caused  which  could  even  require  a  retrograde  maneuver  to  correct  the  semi-major 
axis.  A  simple  way  to  observe  this  phenomenon  in  the  yaw  angle  is  to  examine  the 
velocity  vectors  in  the  normal  and  tangential  planes,  see  Figure  4.7.  In  this  figure,  the 


Figure  4.7.  Velocity  Vectors  in  Normal  and  Tangential  Planes 

velocity  vectors  before  and  after  the  burn,  Vi  and  v2,  must  be  of  equal  magnitude, 
thus  following  the  expression: 


v2  =  Vi  4  Av 

(4.45) 

V2|  =  |vi|=v 

(4.46) 

This  constraint  on  the  velocity  vectors  of  the  satellite  before  and  after  the  burn 
create  an  isosceles  triangle  allowing  the  yaw  angle,  0O,  to  be  determined  with  the  law 
of  cosines: 

,  7T  1  Av  .  . 

■00  =  7T  +  arcsm  - —  (4.47) 

Z  Z  V 

The  placement  of  the  actual  burn  may  also  effect  a  change  in  the  longitude  of  the 

ascending  node,  which  is  not  desired.  This  change  may  be  nulled  by  performing  the 
burn  at  the  equator  crossing,  which  can  be  seen  by  the  equation: 

Aflsinz  =  —  Av^smu  4  0(e2)  (4.48) 

na 
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The  yaw  angle  must  also  account  for  possible  errors.  There  are  factors  such  as  the 
uncertainty  of  the  actual  yaw  angle,  the  yaw  maneuver  bias,  thruster  misalignment, 
attitude  uncertainty,  and  a  deliberate  yaw  bias  to  achieve  a  small  decrease  in  the  semi¬ 
major  axis  with  the  inclination  adjust  maneuver  since  the  inclination  also  affects  the 
nodal  period.  This  bias,  approximately  2°  for  LANDSAT  5,  avoids  an  immediate 
retrograde  maneuver  to  correct  the  semi-major  axis  for  the  new  nodal  period.  The 
total  impact  of  the  uncertainties  affecting  the  yaw  angle,  can  be  determined: 

A  =  ^£(A^)2  (4.49) 

Now  the  actual  yaw  angle,  ip,  can  be  related  to  the  ideal  yaw  angle,  ipo,  with  the 
equation: 

ip  ~  ipo  +  Aip  (4.50) 


Prom  Figure  4.7,  the  magnitude  of  the  impulsive  velocity  needed  to  change  the  incli¬ 
nation  is  determined: 


Av  = 


Avn 
sin  ip 


(4.51) 


This  value  of  the  impulsive  velocity  can  be  used  to  initialize  the  burn  parameters  in 
the  actual  finite  burn  model  for  the  change  in  inclination  [62]. 


Ground  Track  Maintenance 

As  reviewed  previously,  the  ground  track  drift,  in  kilometers,  is  determined  with 
the  following  commensurability  condition  [5]: 

S  =  { ue-  Cl)PNRe  (4.52) 

where  S  is  the  spacing  between  consecutive  ascending  equator  crossings,  u>e  is  the 
mean  rotation  rate  of  the  Earth,  Cl  is  the  rate  of  the  longitude  of  the  ascending  node, 
PN  is  the  nodal  period,  and  Re  is  the  equatorial  radius  of  the  Earth.  According  to 
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McClain  [64],  the  ground  track  drift  after  one  orbit  revolution,  AS,  can  be  displayed 
as  a  truncated  Taylor  series  expansion  if  the  perturbations  axe  assumed  to  be  small; 
this  is  true  if  mean  elements  are  used. 


where 


A  „  SS  A  6S  A  . 
AS  —  —A  a  +  —A  i 
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(4.55) 

(4.56) 

(4.57) 


The  average  rate  of  drift  over  one  day,  also  the  average  rate  of  the  total  ground  track 
drift,  is  determined: 
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(4.60) 
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(4.61) 

(4.62) 


By  integrating  the  expression  for  the  total  ground  track  drift,  the  quadratic  expression 
is  yielded: 

S  —  Sq  =  Ki{t  —  <o)  +  ~d^2{t  —  to)2  (4.63) 
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where  S  and  So  are  the  values  of  the  ground  track  drift  from  the  nominal  ground 
track  at  their  respective  times  of  t  and  t0.  By  allowing  S  to  equal  the  maximum 
longitude,  Smax  in  this  equation,  the  basis  to  calculate  the  necessary  Av  to  cause  the 
ground  track  to  drift  to  its  maximum  longitude,  SmaXi  is  provided.  First  the  partial 
derivatives  of  S  with  respect  to  the  semi-major  axis  and  the  inclination  must  be 
calculated  from  Equations  4.54  and  4.55.  Then  the  value  for  the  quadratic  coefficient 
K2  must  be  calculated  from  Equation  4.62.  From  this  value  of  K2 ,  the  value  for  the 
coefficient  K\  can  be  determined  by  inverting  Equation  4.63: 

K,  =  y/2K2(SQ  -  Smax)  (4.64) 

The  desired  value  for  the  semi-major  axis  adjust  is  determined  by  inverting  Equa¬ 
tion  4.61: 

Ao  =  K1PAr^  (4.65) 

The  impulsive  change  in  velocity  to  change  the  semi-major  axis  will  be  in  the  tan¬ 
gential  direction  and  is  computed: 

2„ 

A  vt  =  — — (4.66) 
2v 

If  the  orbit  eccentricity  is  assumed  to  be  very  small,  as  in  the  LANDSAT  6  case,  the 
impulsive  change  in  velocity  can  be  more  simply  expressed: 

AvT  =  ^Aa  +  0{e)  (4.67) 

This  value  of  the  impulsive  velocity  can  be  used  to  initialize  the  burn  parameters  in 
the  actual  finite  burn  model  for  the  change  in  the  semi-major  axis  [64]. 

Frozen  Orbit  Maintenance 

As  discussed  above,  the  frozen  orbit  condition  is  a  function  of  the  calculated 
nominal  values  of  eccentricity,  cq,  and  argument  of  perigee,  u>o-  An  additional  way  of 
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representing  these  parameters  is  to  combine  them  to  represent  an  eccentricity  vector 
in  polar  coordinates,  where  [61]: 

ex  =  ecosu;  (4.68) 

ey  =  esinu>  (4.69) 


The  actual  eccentricity  vector  will  be  displaced  from  the  nominal  eccentricity  vector 
because  of  orbit  perturbations  that  were  not  modeled  in  the  frozen  orbit  condition. 
The  actual  eccentricity  vector  will  be  displaced  from  the  nominal  eccentricity  vector 
by  the  vector  Ae  with  components  [63]: 


Ae*  =  e  cos  u  —  e0  cos 
Aey  =  esincj  —  eosinc^o 

The  magnitude  of  the  displacement  eccentricity  vector  is  calculated: 

cP  =  |Ae| 

=  \j  (Aez)2  +  (Ae„)2 

and  the  direction  of  the  displacement  eccentricity  vector  is  calculated: 

'  Ae„ 


(4.70) 

(4.71) 


(4.72) 

(4.73) 


0  =  arctan 


VAeJ 


(4.74) 


These  polar  coordinates  of  the  eccentricity  vector  are  sensitive  to  changes  in  the 
velocity  in  the  tangential  direction.  The  effects  of  these  changes  are  determined  [61]: 


6ex  =  — cosuA  vt  +  O(e) 
na 

2 

6ev  =  — sinuAur  +  O(e) 
na 


(4.75) 

(4.76) 


where  n  is  the  satellite  mean  motion,  a  is  the  semi-major  axis,  and  u  is  the  argument 
of  latitude  equal  to  the  sum  of  the  true  anomaly  and  the  argument  of  perigee.  The 
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argument  of  latitude  where  the  burn  takes  place  determines  the  effect  that  the  burn 
will  have  on  the  eccentricity  vector.  As  observed  from  Equation  4.67,  the  tangential 
velocity  needed  to  adjust  the  semi-major  axis  is  independent  of  the  latitude  of  the 
burn.  Thus,  the  burn  may  be  planned  in  its  location  to  adjust  not  only  the  semi¬ 
major  axis  but  also  the  frozen  orbit  parameters  back  to  their  nominal  values  with 
no  more  fuel  expended  than  was  previously  necessary  by  the  independent  semi-major 
axis  adjust. 

The  magnitude  of  the  displacement  of  the  eccentricity  vector  can  be  combined 
from  Equations  4.75  and  4.76  in  the  following  expression  [63]: 

D2  =  |<5e|2  (4.77) 

=  (— )\avt?  (4.78) 

\naj 

where  the  value  of  the  change  in  velocity  in  the  tangential  direction,  A vt  is  determined 
from  the  needed  change  in  velocity  for  the  semi-major  axis  adjust.  If  the  needed 
change  in  the  eccentricity  vector,  d,  is  less  than  the  change  in  the  eccentricity  vector 
due  to  the  change  in  tangential  velocity,  D,  then  two  separate  burns  are  needed  to 
correctly  adjust  both  the  semi-major  axis  and  the  frozen  orbit  parameters.  This  is 
called  Case  1.  If  the  needed  change  is  greater  than  or  equal  to  the  change  in  the 

eccentricity  vector  due  to  the  change  in  tangential  velocity,  only  one  orbit  burn  is 

needed.  This  is  called  Case  2  [63].  In  Case  1,  the  two  burns  will  take  place  at: 

Ui  =  6  ±  7r  (4-79) 

«2  =  lii  +  7T  (4.80) 

where  u  is  the  argument  of  latitude  where  the  burn  is  to  take  place  and  0  is  defined 
in  Equation  4.74.  The  magnitude  of  impulsive  velocity  expended  at  each  of  the  two 
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burns  is  determined: 


A„Ti  _  ^L  +  ^d  (4,81) 

**  -  («2) 

where  d  is  determined  by  Equation  4.73,  A vt  is  determined  from  Equation  4.67,  n  is 
the  satellite  mean  motion,  and  a  is  the  semi-major  axis.  In  Case  2,  the  one  impulsive 
burn  will  take  place  at  an  argument  of  latitude: 

u  =  9  +  7T  (4.83) 

where  0  is  determined  in  Equation  4.74.  The  magnitude  of  the  burn  in  Case  2  is 
equal  to  the  impulsive  burn  calculated  in  Equation  4.67  for  the  change  in  the  semi¬ 
major  axis.  These  values  of  the  impulsive  velocities  can  be  used  to  initialize  the  burn 
parameters  in  the  actual  finite  burn  model  for  the  change  in  the  semi-major  axis  [63] 

4.2.2  Finite  Burn  Models 


As  discussed  in  the  previous  section,  the  LANDSAT  6  mission  profile  has  finite 


burn  models  that  are  to  be  initialized  by  the  impulsive  changes  in  velocity  to  adjust 
the  orbital  elements  to  maintain  the  repeat  ground  track,  sun  synchronous,  and  frozen 


orbit  conditions.  The  finite  burn  model  over  time  is  expressed  [59]: 


jg*:(i  -ig)  _  j,  j 

A{\n[L{t  -  to)]  +  C) 
1 

e-d(t-tb) 


0  <  t  -  ty  <  40  msec 
40  msec  <  t  —  to  <  867  msec 
867  msec  <  t  —  t0  <  tb  —  t0 
0  <  t  —  tb  <  104  msec 


where 


(4.84) 


e(t )  =  thrust  acceleration  magnitude  normalized  by  the  gravitational 
acceleration  at  the  Earth’s  surface 


e0  =  maximum  value  of  the  normalized  thrust  acceleration 
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to  =  burn  acceleration  initiation  epoch 


tb  =  burn  acceleration  cut-off  epoch 

k,  A,  L,C,d  =  thruster  dependent  parameters 

This  model  for  the  thrust  acceleration  may  be  added  as  a  perturbation  to  the  equa¬ 
tions  of  motion  and  integrated  either  analytically  or  numerically  to  determine  how 
the  mean  elements  change  with  respect  to  time  during  the  maneuver.  Numerical 
integration  would  require  step  sizes  on  the  order  of  fractions  of  a  second.  This  is 
inconsistant  with  the  general  approach  in  the  semianalytical  theory.  An  analytical 
integration  would  assume  an  upper  limit  of  16  minutes  to  the  burn  time,  or  ap¬ 
proximately  60°  of  orbit  motion.  These  analytical  models  were  applied  to  a  finite 
burn  lasting  longer  than  the  16  minute  limit,  ie.  18  min.  and  30  min.  These  burns 
calculated  a  semi-major  axis  that  differed  from  the  more  accurate  numerical  integra¬ 
tion  scheme  by  2.5  x  10-3%  for  the  18  minute  burn  and  6  x  10-3%  for  the  thirty 
minute  burn.  These  errors  are  very  small  [66].  The  analytical  equations  used  in  the 
LANDSAT  6  mission  design  for  the  actual  thrust  profile  are  determined  by  McClain 
in  references  [67,59]. 
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Chapter  5 

Ground  Track  Motion  Forecasting 


Forecasting  the  ground  track  is  a  difficult  task.  An  orbit  propagator  can  be  used 
to  determine  to  a  certain  accuracy  the  perturbations  affecting  the  orbit  and  how  these 
changes  affect  the  ground  track.  The  most  difficult  of  the  perturbations  to  forecast 
is  the  atmospheric  drag.  Complex  equations  for  drag  dependent  on  the  spacecraft 
model  exist,  but  a  simple  equation  for  determining  atmospheric  drag  is  [23]: 

_  1  Ct>A  _  _  /E  lN 

&drag  —  2  TYl  ^re^re^  (5-1) 

where  adrag  is  the  perturbing  acceleration  due  to  drag,  Co  is  the  coefficient  of  drag 
dependent  on  the  spacecraft  structure  and  material,  A  is  the  area  of  the  satellite 
perpendicular  to  the  satellite  velocity  relative  to  the  atmosphere,  vre/,  and  m  is 
the  spacecraft  mass.  The  upper  atmospheric  density,  p ,  is  affected  by  a  complex 
interaction  between  the  sun,  the  solar  wind,  and  the  Earth’s  geomagnetic  field.  Thus, 
the  atmospheric  density  is  a  function  on  the  satellite’s  altitude,  latitude,  longitude, 
local  time,  as  well  as  geomagnetic  and  auroral  activity,  and  solar  flux  levels  [75].  The 
solar  flux  affects  the  atmospheric  density  through  direct  and  instantaneous  heating 
by  extreme  ultraviolet  radiation.  The  geomagnetic  activity  affects  the  atmosphere 
through  delayed  indirect  heating  of  atmospheric  energetic  particles  from  collisions 
with  charged  particles  emitted  from  the  sun.  This  heating  of  the  atmosphere  causes 
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the  density  to  increase  at  higher  altitudes.  The  solar  flux  level  and  the  geomagnetic 
activity  are  difficult  to  predict  but  are  two  main  inputs  into  current  models  of  the 
atmosphere.  These  two  factors  are  discussed  further  below. 

5.1  Solar  Activity 

As  stated  in  Prochaska  [75],  solax  flux,  the  amount  of  heating  that  the  upper 
atmosphere  experiences  through  the  absorption  of  extreme  ultraviolet  (euv)  solar 
radiation,  is  impossible  to  measure  since  the  atmosphere  does  not  allow  the  euv 
radiation  to  pass  through.  There  are  no  current  space  born  systems  to  measure 
euv  flux  nor  atmospheric  density  models  to  use  these  space  based  measurements. 
Investigation  of  new  models  and  new  parameters  such  as  the  precipitation  index 
is  currently  being  done  [37,38].  Therefore,  the  euv  measurement  must  be  inferred 
from  the  closely  correlated  Earth  based  measurements  of  10.7  cm  length  radio  waves 
(Fl0.7).  Both  radiation  types,  euv  and  F10.7,  originate  in  the  same  layers  of  the  sun’s 
chromosphere.  The  Earth’s  atmosphere  is  transparent  to  the  F10.7  radiation;  its  value 
has  been  routinely  measured  since  1940.  The  F10.7  radiation  is  measured  in  solar  flux 
units  (SFU),  1  SFU  =  1  x  10~22 ;  typical  values  range  from  less  than  70  to 
more  than  300  SFU  [74]. 

The  solar  flux  exhibits  two  superimposed  cyclic  variations.  The  primary  factor  is 
an  approximate  11  year  cycle  that  roughly  parallels,  but  lags  a  few  years  behind  the 
sun  spot  cycle.  The  minimum  of  this  cycle  is  not  halfway  between  the  two  maximums 
since  the  decreasing  phase  of  the  cycle  is  6-7  years.  The  actual  peak  of  the  11  year 
cycle  varies  from  cycle  to  cycle. 

The  secondary  solar  flux  period  lasts  approximately  6  months  and  is  related  to 
the  varying  distance  of  the  Earth  from  the  sun  during  its  slightly  elliptical  orbit.  In 
addition  to  these  two  cycles,  there  are  irregular  changes  to  the  solar  flux  that  are 
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related  to  the  growth  and  decay  of  active  solar  regions.  These  active  regions  not  only 
have  many  different  patterns  of  growth,  stability,  and  decay,  but  also  rotate  with  the 
sun,  approximately  27  days  per  rotation,  thus  causing  an  uncertain  cyclic  pattern. 
These  patterns  are  difficult  to  predict.  It  is  also  difficult  to  determine  how  they  will 
contribute  to  the  heating  of  the  Earth’s  atmosphere.  It  is  this  uncertainty  that  affects 
the  accuracy  of  solar  flux  predictions  [74], 

Another  characteristic  of  the  density  is  the  diurnal  bulge,  or  local  atmospheric 
expansion,  created  by  the  constant  heating  of  the  atmosphere  on  the  sunlit  side  of 
the  Earth.  The  bulge  axis  has  a  local  time  of  approximately  2-2:30  pm.  The  bulge 
is  centered  on  the  equator  at  the  equinoxes  but  moves  to  higher  latitudes  than  the 
sun’s  declination  at  the  solstices.  The  diurnal  bulge  makes  the  atmospheric  density 
dependent  on  latitude,  local  time,  and  time  of  year  [23] 

The  most  commonly  accepted  measurement  of  F10.7  is  distributed  by  the  National 
Oceanic  and  Atmospheric  Administration  (NOAA)  at  the  National  Geophysical  Data 
Center  in  Boulder,  Colorado.  The  actual  measurement  is  made  at  the  Algonquin 
Radio  Observatory  in  Ottawa,  Ontario,  Canada  [82]. 

5.2  Geomagnetic  Index 

Direct  collisions  of  the  solar  wind  with  air  particles  interacting  with  the  Earth’s 
geomagnetic  field  heat  the  atmosphere.  Geomagnetic  activity  must  be  measured  to 
determine  the  heat  generated.  The  planetary  geomagnetic  index,  kp,  is  a  quasi- 
logarithmic  worldwide  average  of  geomagnetic  activity  below  the  auroral  zones.  This 
value  is  measured  every  three  hours.  The  geomagnetic  planetary  amplitude,  op,  is 
the  linear  equivalent  of  this  index.  Eight  values  are  measured  daily  and  averaged  to 
create  the  commonly  used  daily  planetary  amplitude,  Ap  [23,36].  The  daily  planetary 
amplitude  Ap  is  measured  in  units  of  2  gammas,  where  1  gamma  =  10-9  Tesla.  A 
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typical  range  of  Ap  is  from  0  to  400;  values  greater  than  100  are  rare.  The  daily 
planetary  amplitude  tends  to  follow  the  11  year  sun  spot  cycle,  although  there  are 
consistently  large  maximums  of  Ap  in  the  declining  phase  of  each  11  year  cycle.  There 
is  also  a  secondary  semi-annual  cycle  due  to  the  variable  position  of  the  solar  wind 
with  respect  to  the  Earth’s  magnetosphere.  This  cycle  is  just  as  variable  and  as  hard 
to  predict  as  the  sun  spot  cycle.  Variations  of  Ap  from  the  sun  spot  and  semi-annual 
cycles  are  mainly  due  to  solar  flares,  coronal  holes,  disappearing  solar  filaments,  and 
the  near  Earth  solar  wind  environment  [36,74].  Intense  geomagnetic  activity  at  the 
auroral  zones  affects  the  shape  of  the  atmosphere  and  makes  atmospheric  density 
dependent  on  latitude  [23]. 

Planetary  geomagnetic  indices  kp  and  ap  are  compiled  using  measurements  from 
eleven  observatories  which  lie  between  46  N  and  63  S  latitudes;  3  of  these  in  the 
United  Kingdom,  2  are  in  Canada,  2  are  in  the  USA,  and  the  remaining  four  in  New 
Zealand,  Australia,  Sweden,  and  Denmark.  The  most  accepted  compilation  of  the 
measurements  from  these  observatories  is  created  by  the  Institut  fur  Geophysik,  at 
Gottingen  University,  Germany  [39]. 

5.3  Confidence  in  Forecast 

Upper  and  lower  confidence  limits  for  the  acceleration  due  to  atmospheric  drag 
are  required  to  determine  a  conservative  time  and  magnitude  of  an  orbit  maintenance 
burn.  When  calculating  the  epoch  of  a  maintenance  burn,  it  is  good  to  assume  an 
upper  bound  for  drag  to  ensure  that  the  ground  track  of  the  satellite  does  not  go 
beyond  the  lower  ground  track  boundary  before  the  maintenance  burn  takes  place. 
In  addition,  it  is  wise  to  use  a  lower  bound  of  the  atmospheric  drag  to  determine  the 
magnitude  of  the  burn.  If  the  actual  atmospheric  density  is  less  than  the  density  used 
in  calculating  the  magnitude  of  the  burn,  then  the  ground  track  will  exceed  its  upper 
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boundary  and  an  unplanned  and  unwanted  retrograde  burn  would  be  required  to  set 
the  ground  track  within  its  nominal  bounds. 

As  seen  in  Equation  5.1,  the  acceleration  due  to  drag  is  proportional  to  the  co¬ 
efficient  of  drag,  Cd,  and  the  atmospheric  density,  p.  In  principle,  confidence  limits 
of  the  forecasted  atmospheric  drag  could  be  inferred  from  statistics  on  Fw  7  and 
Ap,  but  this  approach  is  difficult  since  atmospheric  models  axe  complicated  and  the 
independence  of  the  variables  F10.7  and  Ap  can  not  be  assumed.  Research  for  the 
TOPEX/POSEIDON  mission  assumes  that  kp  variations  are  independent  of  *10.7 
variations.  Reference  [6]  infers  this  from  observing  actual  values  for  the  solar  flux 
and  geomagnetic  activity.  It  is  not  stated,  however,  if  this  assumption  is  possible 
because  of  the  high  altitude  of  TOPEX.  In  an  alternative  method,  nominal  density 
predictions  could  be  used  together  with  upper  and  lower  estimates  for  the  drag  co¬ 
efficient;  Co  could  be  varied  using  1  o  estimates  from  current  orbit  determination 
or  longer  term  statistics  on  Co-  Neither  of  these  methods  seems  very  reliable.  The 
method  used  in  this  thesis  to  determine  upper  and  lower  confidence  limits  of  atmo¬ 
spheric  drag  is  to  use  long  term  statistics  of  the  atmospheric  density  as  modeled  to 
determine  upper  and  lower  estimated  limits  of  the  forecasted  density. 

5.3.1  Soviet  Density  Model 

The  LANDSAT  6  mission  will  use  a  Jacchia-Roberts  density  model  in  determining 
the  maintenance  burns  needed  for  stationkeeping.  Currently,  there  is  not  a  standalone 
program  for  Jacchia-Roberts  density  evaluation  at  Draper  Laboratory.  A  simple 
Soviet  density  model  [85]  is  used  in  this  thesis  to  determine  a  first  approximation 
to  the  density  confidence  interval.  This  Soviet  density  model  was  successfully  used 
in  previous  LANDSAT  6  research.  The  variability  in  the  Soviet  model  gave  a  ler 
prediction  error  similar  to  the  position  error  of  actual  data  calculated  with  a  IIarris- 
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Priester  density  model  [13].  This  Soviet  model  was  also  compared  to  the  Jacchia  1971 
density  model  and  found  to  differ  only  5-10%  [86].  A  future  study  will  employ  the 
method  used  here  with  the  Jacchia-Roberts  density  model.  The  difference  between 
the  confidence  interval  determined  here  and  that  determined  by  the  Jacchia-Roberts 
model  is  not  anticipated  to  exceed  10%. 

The  Soviet  model  was  constructed  from  observations  of  the  orbit  motion  of  the 
Soviet  Kosmos  satellites.  This  model  includes  the  dependence  of  the  density  on 
solar  flux  and  geomagnetic  activity  as  well  as  the  diurnal  and  semi-annual  density 
variations.  This  model  is  recommended  for  satellites  in  the  altitude  range  of  160-600 
km.  Since  the  LANDSAT  6  mission  altitude  is  700  km,  the  application  to  this  mission 
is  necessarily  tentative. 

The  Soviet  formula  models  the  actual  atmospheric  density  as  a  product  of  five 
factors,  where  each  factor  corresponds  to  a  particular  density  variation  [85]: 

p  =  pnk\k2k^k4  (5.2) 

where  pn  is  the  night  time  density  profile  assumed  to  be  exponentially  decreasing 
with  respect  to  altitude: 

pn  =  exp  [a  i  -  a2[h  -  a3)*]  (5.3) 

The  coefficients  in  this  expression  were  determined  empirically  and  are  listed  in  Ta¬ 
ble  5.1.  The  altitude  of  the  satellite  in  km  is  listed  as  the  variable  h. 

The  k\  factor  in  Equation  5.2  represents  the  dependence  of  the  atmospheric  density 
on  the  solar  flux  *1 0.7-  The  factor  k2  accounts  for  the  diurnal  effect  and  k3  takes  into 
account  the  semi-annual  effect  of  the  density  variation.  The  factor  k4  represents  the 
effect  of  geomagnetic  activity  index  ap  [85].  In  determining  the  confidence  interval, 
only  daily  averages  of  the  density  are  needed  and  the  dependence  on  the  semi-annual 
cycle  is  not  needed  if  the  error  of  the  predicted  density  is  normalized  by  the  predicted 
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nominal  density.  This  normalization  not  only  makes  the  density  a  stationary  process, 
but  also  removes  the  need  for  the  inclusion  of  the  night-time  density,  pn ,  in  the 
calculations  [13].  The  k t  and  k 4  factors  given  below  were  the  only  factors  used  in  this 
thesis: 

,  ,  ,  (&i  +  b2h)(Fi 0.7  —  F) 

kl  =  i  + - ^ - 

k±  =  1  +  (ei  +  e2h)  In 

where  the  e  and  6  coefficients  are  listed  Table  5.1.  The  variable  F  represents  the 
mean  solar  flux  for  the  period  looked  at.  The  Soviet  density  model  only  allows  for  a 
mean  solar  flux  of  75,  100,  125,  or  150.  Thus,  the  actual  average  solar  flux  must  be 
rounded  to  one  of  these  values.  The  term  ap  represents  the  mean  geomagnetic  index 
dependent  on  the  mean  solar  flux  which  is  shown  in  Table  5.1. 


(3.4) 

(5.5) 


Table  5.1.  Soviet  Density  Model  Parameters 


■ 

Value  with  F  = 

75 

100 

125 

150 

a  1 

-14.030 

-15.095 

-17.028 

-16.072 

a2 

0.9108 

0.8299 

0.7198 

0.7155 

03 

59.77 

68.92 

93.36 

70.33 

'1® 

-0.630 

-0.710 

-0.765 

62 

0.00506 

0.00560 

0.00562 

0.00571 

ei 

-0.132 

-0.130 

-0.128 

-0.115 

e2 

0.00108 

0.00104 

0.00095 

0.00089 

do 

2 

2 

3 

4 
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5.3.2  Algorithm  for  Upper  and  Lower  Bounds  of  Interval 


In  the  Draper  version  of  the  Goddard  Trajectory  Determination  System  (Draper 
GTDS),  the  magnitude  of  the  acceleration  due  to  drag  is  calculated: 

1  A 

CLdrag  =  »CD(  1  +  t)p—  \vrel  |2  (5.6) 

z  m 

where  the  variables  are  the  same  as  in  Equation  5.1  with  the  exception  of  e  which  is 
a  correction  factor.  Often  orbit  determination  runs  are  done  with  the  coefficient  of 
drag,  Cq,  held  constant  and  the  density,  p,  given  by  a  fixed  atmospheric  model;  the 
correction  factor,  e,  which  is  nominally  zero,  can  be  solved  for  in  the  differential  cor¬ 
rection  to  better  approximate  the  actual  density  encountered.  The  proposed  mission 
planning  software  for  LANDS  AT  6  will  use  the  variable  e  as  an  instrument  to  adjust 
the  density  to  the  upper  and  lower  confidence  limits. 


Pup  (1  T  Cup) Pnom 

(5.7) 

Plow  (1  T  Clow) Pnom 

(5.8) 

where  pn0m  is  the  nominal  predicted  density,  tup  and  t[OW  are  the  upper  and  lower 
confidence  limits  for  the  normalized  density  error ,  and  pup  and  piow  are  the  resulting 
upper  and  lower  confidence  limits  for  the  density.  Note  that  is  is  assumed  that  the 
value  e;ou,  is  negative. 

The  normalized  density  error  is  calculated: 


Pact  Pi 
Pi 


(5.9) 


where  pi  is  the  atmospheric  density  calculated  from  the  predicted  solar  flux  and 
geomagnetic  activity  indices  for  a  particular  number  of  days  in  the  future,  /,  and  pact 
is  the  density  calculated  from  the  actual  measured  solar  flux  and  geomagnetic  data  for 
that  day.  The  National  Oceanic  and  Atmospheric  Administration  (NOAA)  publishes 
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the  predicted  solar  flux  and  geomagnetic  activity  index  pairs  weekly  for  up  to  27  days 
in  the  future.  Thus,  the  possible  latency  of  a  forecast,  l,  ranges  from  1  to  27.  The 
predicted  values  of  the  solar  flux  can  often  be  normalized  to  account  for  the  varying 
distance  of  the  Earth  from  the  sun.  This  thesis  uses  un-normalized  values  of  the  solar 
flux  since  the  omission  of  the  k3  factor  in  the  Soviet  density  model  averages  out  that 
semi-annual  cycle  and  it  was  decided  that  double  averaging  would  not  ameliorate  the 
calculated  output.  It  can  be  seen  from  Equation  5.9  that  the  absolute  lower  limit  of 
the  normalized  error  is  negative  one  and  the  absolute  upper  limit  is  not  defined. 

To  determine  the  density  confidence  interval,  the  upper  bound  would  be  set  so 
that  some  percentage,  P,  of  the  normalized  errors  lie  below  that  bound.  The  lower 
bound  would  be  set  so  that  the  same  percentage,  P,  would  lie  above  that  bound. 
This  method  does  not  ensure  that  the  bounds  will  be  symmetric  about  the  mean 
normalized  error.  In  LANDSAT  6,  the  percentage,  P,  is  required  to  be  97.3%  and  is 
the  value  used  here  in  this  thesis. 

One  undetermined  variable  of  the  method  discussed  above  is  the  length  of  time 
over  which  the  density  will  be  looked  at  to  determine  the  upper  and  lower  bounds. 
A  long  sample  period  is  recommended  to  increase  the  reliability  of  the  statistics,  but 
the  solar  cycle  upon  which  the  solar  flux  and  geomagnetic  indices  are  dependent  is 
not  a  stationary  process,  ie.  the  error  of  the  predictions  will  change  over  the  time  of 
the  cycle.  In  this  non-stationary  case,  the  use  of  a  small  sample  size  is  recommended. 

Method  1,  which  assumes  that  the  normalized  error  creates  a  stationary  solar  cy¬ 
cle,  examined  solar  flux  and  geomagnetic  data  from  1986-1990.  For  the  non-stationary 
case,  Method  2,  only  data  for  1990,  the  most  recent  completed  previous  year,  was 
analyzed,  since  this  year  has  the  closest  ties  to  1992,  the  first  year  of  the  LANDSAT 
6  mission. 
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5.4  Results 


The  actual  solar  flux  and  all  predicted  values  were  received  from  the  Jet  Propulsion 
Laboratory  (JPL).  The  actual  solar  flux  values  were  compiled  for  the  years  1986-1990 
from  the  Preliminary  Report  and  Forecast  of  Solar  Geophysical  Data  (PRF)  published 
weekly  by  the  National  Oceanic  and  Atmospheric  Administration  (NOAA).  Data 
was  given  for  previous  years  but  the  year  1985  was  missing  many  data  points  due 
to  magnetometer  outages,  therefore  1984  and  1985  were  not  included.  The  27  day 
forecast  solar  flux  data  is  also  predicted  weekly  by  NOAA  in  the  same  publication 
and  was  compiled  by  JPL.  The  actual  values  for  the  geomagnetic  daily  planetary 
amplitude,  Ap,  used  in  this  thesis  were  received  from  the  National  Geophysical  Data 
Center,  a  division  of  NOAA.  This  data  was  deemed  more  accurate  than  the  data  from 
JPL,  since  the  Ap  averages  from  NOAA  included  measurements  from  remote  stations 
and  the  geomagnetic  index  used,  Ap,  was  consistent  throughout  all  of  the  years.  The 
information  from  JPL  began  with  actual  values  for  the  A  index,  which  is  determined 
from  the  high  latitude  observations  from  Anchorage.  This  index  is  similar  but  not 
identical  to  the  planetary  amplitude,  Ap,  needed  by  the  Soviet  density  model.  This 
data  was  placed  in  files  as  flat  matrices,  ie.  the  same  number  of  columns  in  each  row, 
for  access  by  a  MATLAB  script  file. 

The  Soviet  density  model  was  programmed  as  a  function  to  be  used  by  a  MATLAB 
scripc  file. 

5.4.1  Method  1 

This  method  assumes  that  the  normalized  error  is  stationary  over  the  solar  cycle. 
Therefore,  more  reliable  statistics  can  be  taken  by  using  a  longer  time  period.  In  this 
thesis,  all  of  the  available  data  was  used,  ie.  1986  -  1990.  This  method  determines 
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the  normalized  errors  for  all  of  the  data  and  then  chooses  the  upper  and  lower  bounds 
for  each  latency  day  to  match  the  required  percentages. 

There  was  a  small  problem  with  the  Soviet  density  model  in  the  year  of  1988.  This 
year  started  with  fairly  low  values  of  the  solar  flux  and  then  reached  extremely  high 
values  towards  the  end  of  the  year  as  the  solar  maximum  approached.  This  caused 
negative  values  of  density  to  be  calculated  in  the  beginning  of  1988.  Hence,  1988  was 
split  and  the  densities  calculated  for  each  half  using  the  appropriate  mean  solar  flux 
value  for  that  half.  This  resolved  the  problem. 

The  graph  of  the  normalized  error  limits  can  be  viewed  in  Figure  5.1.  This  data 
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Figure  5.1.  1986-1990  Upper  and  Lower  Limits  of  Confidence  Interval 


is  very  noisy  but  the  general  trend  increases  as  the  latency  of  the  forecast  increases. 
This  general  trend  can  be  viewed  in  Figure  5.3. 
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5.4.2  Method  2 


The  second  method  uses  the  most  recent  year  of  available  solar  flux  and  geo¬ 
magnetic  data,  1990,  to  determine  the  difference  between  the  predicted  and  actual 
densities.  These  errors  are  normalized  and  then  the  confidence  interval  determined. 
These  bounds  could  then  be  used  in  Equations  5.7  and  5.8. 

This  method  assumes  that  for  missions  like  LANDSAT  6,  where  the  11  year  solar 
cycle  is  declining,  the  previous  year  of  data  would  provide  a  conservative  estimate  of 
the  current  density  confidence  interval.  This  method  also  assumes  that  the  normalized 
error  is  accurate  for  only  one  year,  which  is  more  accurate  than  the  assumption  in 
Method  1.  Method  2  should  be  repeated  approximately  yearly  throughout  the  mission 
lifetime  to  determine  the  new  confidence  interval  limits  using  data  from  the  previously 
completed  year  so  that  the  method  assumptions  are  not  violated  and  the  best  results 
achieved.  In  actuality,  this  could  be  repeated  at  anytime  when  it  is  determined  that 
new  limits  are  needed. 

The  upper  and  lower  limits  of  the  confidence  interval  can  be  viewed  in  Figure  5.2. 
This  data  is  even  more  noisy  than  that  calculated  for  1986-1990  in  Method  1.  This 
was  to  be  expected  since  1990  is  near  the  solar  maximum  and  the  errors  in  the  forecast 
are  more  extreme.  Again,  the  general  trend  of  the  confidence  interval  is  to  increase 
as  the  latency  of  the  forecasted  data  increases  as  was  expected.  This  can  be  viewed 
in  Figure  5.3. 

5.5  Conclusions 

As  stated  previously,  the  results  of  both  methods  can  be  compared  in  Figure  5.3. 

The  actual  equations  for  these  general  trends  can  be  found  in  Table  5.2.  where 
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Figure  5.2.  1990  Upper  and  Lower  Limits  of  Confidence  Interval 
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Figure  5.3.  Upper  and  Lower  Limits  of  Confidence  Interval 
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Table  5.2.  Equations  for  Trends  of  Upper  and  Lower  Limits 


1986-1990 

eup  =  0.77288  +  2.3144  x  10 ~'H 
ei^  =  -0.62798  -  4.5675  x  10~3/ 

1990 

eup  =  0.78666  +  8.9609  x  10~3f 
etow  =  -0.66662  -  5.5858  x  10~3Z 

l  is  the  latency  value  of  the  forecast,  1  <  l  <  27.  These  general  trends  are  better 
estimations  of  the  actual  upper  and  lower  limits  of  the  density  since  they  axe  smoother 
and  are  not  dependent  on  the  exact  error  that  might  be  specific  to  a  certain  year. 

The  confidence  interval  found  for  the  cumulative  years  of  1986-1990  in  Method  1 
gives  a  slightly  more  conservative  est  mate.  Although  the  lower  limit  for  Method  1,  is 
slighlty  higher  than  the  lower  limit  for  Method  2,  their  difference  is  negligible  and  the 
lower  limits  for  both  methods  are  approximately  equal.  Therefore,  it  is  recommended 
that  the  1986-1990  limits  define  the  confidence  interval  implemented  for  LANDSAT 
6.  If  it  is  found  during  the  mission  lifetime  that  these  limits  are  too  conservative,  the 
results  from  Method  2  could  be  substituted.  It  is  better  to  staxt  out  the  mission  with 
a  limit  that  is  too  conservative  requiring  many  orbit  adjust  burns  and  then  switch  to 
a  looser  estimate.  If  the  initial  confidence  interval  used  is  too  narrow,  the  bounds  of 
the  ground  track  could  be  violated  and  the  mission  endangered. 

Because  of  the  makeup  of  the  Draper  Semianalytic  Theory  Standalone  Propagator, 
it  is  desired  that  the  upper  and  lower  limits  be  a  constant  rather  than  varying  with 
the  forecast  latency.  In  this  case,  it  is  advised  that  the  limits  for  the  27th  day  of  the 
forecast  be  used,  since  these  values  are  the  most  conservative.  Again,  these  may  be 
changed  if  they  axe  found  to  be  too  conservative  during  the  mission  lifetime. 

One  foreseeable  problem  with  choosing  the  most  conservative  confidence  interval 
is  that  the  lower  limit  might  equal  approximately  zero,  causing  the  lower  limit  of  the 
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density  to  be  no  atmospheric  density.  This  could  result  in  convergence  difficulties 
in  the  Draper  Semianalytic  Theory  Standalone.  If  this  occurs,  it  is  suggested  that 
a  slightly  increased  lower  limit  be  selected.  But  in  the  analyses  above,  this  problem 
was  not  encountered,  thus  the  actual  lower  limit  is  recommended. 
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Chapter  6 

Results  and  Conclusions 


This  thesis  provides  a  starting  point  for  research  and  development  involving  ma¬ 
neuver  planning  with  a  semianalytic  satellite  theory.  The  common  orbit  constraints  of 
sun  synchronous,  repeat  ground  track,  and  frozen  orbit  were  discussed  in  the  context 
of  maneuver  planning.  Since  drag  is  the  most  uncertain  acceleration  that  perturbs 
these  constraints  from  their  nominal  position,  >.  confidence  interval  for  the  atmo¬ 
spheric  density  was  determined  u.  :ng  solar  flux  and  geomagnetic  data  in  a  simple  but 
fairly  accurate  density  mode’.  This  chapter  concludes  this  thesis  by  recommending  a 
density  confidence  interval  for  use  in  a  semianalytic  satellite  theory  orbit  propagator. 
In  addition,  several  topics  are  suggested  for  further  research. 

6.1  Maneuver  Planning  Software  Tools 

As  discussed  in  Chapter  2,  a  semianalytic  propagation  theory  combines  the 
speed  of  an  analytic  theory  but  retains  the  accuracy  of  a  numerical  theory.  Thus, 
the  semianalytic  satellite  theory  serves  as  an  excellent  candidate  for  application  to 
a  maneuver  planning  process.  The  Draper  Semianalytic  Satellite  Theory  Standalone 
Propagator  is  not  only  fast  and  accurate,  but  it  is  portable  to  many  computing 
environments. 
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Program  MEANELT  is  an  excellent  tool  for  maneuver  planning  but  has  a  few 
problems  if  the  satellite  is  not  in  a  geosynchronous  orbit.  This  difficulty  occurs  if  the 
maintained  longitude  of  the  ascending  node  needs  to  be  calculated  more  than  once 
per  day,  as  in  a  low  altitude  satellite.  This  problem  for  repeat  ground  track  orbits 
may  be  avoided  by  maintaining  the  semi-major  axis  which  is  the  primary  variable  in 
the  longitude  of  the  ascending  node  rate. 

It  is  recommended  that  analytic  equations  be  placed  in  the  Draper  Semianalytic 
Satellite  Theory  to  model  the  variable  thrust  throughout  the  entire  burn.  These 
models  will  be  initialized  with  an  impulsive  burn  and  then  optimized  so  that  the 
actual  burn  achieves  the  desired  change  in  the  orbital  elements.  Expressions  for 
these  impulsive  burn  parameters  can  be  found  in  Chapter  4  and  the  analytic  thrust 
equations  can  be  found  in  references  [67,59]. 

6.2  Density  Confidence  Interval 

In  forecasting  the  ground  track  drift,  the  most  uncertain  perturbing  acceleration 
is  drag.  Therefore,  in  planning  an  orbit  adjust  burn,  the  upper  limit  of  drag  must 
be  used  to  plan  when  to  start  the  burn  so  that  the  eastern  ground  track  limit  is 
not  exceeded.  The  lower  limit  for  drag  must  be  used  to  determine  the  magnitude  of 
the  burn  so  that  the  western  ground  track  limit  is  not  exceeded.  It  is  recommended 
that  atmospheric  drag  be  adjusted  using  upper  and  lower  limits  of  the  atmospheric 
density  determined  from  the  calculated  normalized  error  from  all  the  available  data 
from  previous  years.  In  this  case,  the  years  1986-1990  were  used.  The  upper  and 
lower  limits  will  change  for  the  number  of  days  in  the  future  that  the  predicted  values 
of  solar  flux  and  geomagnetic  activity  forecast;  this  is  called  the  forecast  latency.  This 
change  in  the  limits  for  each  latency  can  be  fairly  noisy,  see  Figures  5.1  and  5.2.  The 
general  trend  of  the  limits  for  the  year  of  1986-1990  can  be  observed  in  Figure  5.3.  The 
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equations  for  this  general  trend  are  listed  in  Table  5.2.  Since  the  Draper  Semianalytic 
Standalone  Propagator  requires  a  constant  value,  the  most  conservative  estimate  is 
desired.  Thus,  the  27  day  latency  limits  for  1986-1990  axe  recommended.  If  these 
values  are  found  to  be  too  conservative  in  estimating  the  atmospheric  drag  on  the 
satellite  during  the  mission,  they  may  be  lowered  to  the  limits  determined  for  the  year 
of  1990  alone.  It  is  better  to  start  out  too  conservative  and  adjust  to  lower  values  to 
ensure  that  the  ground  track  boundaries  are  not  exceeded. 

6.3  Suggested  Further  Research 

The  accuracy  and  successful  application  of  the  Draper  Semianalytic  Standalone 
Propagator  has  been  established.  But  now,  it  must  be  applied  to  the  actual  LANDSAT 
6  mission  to  prove  its  success  at  real  time  maneuver  planning.  In  addition,  the  ap¬ 
plication  of  repeat  ground  track  and  frozen  orbits  to  multiple  satellites  should  be 
investigated  as  a  simple  and  accurate  means  of  constellation  management. 

Program  MEANELT  needs  extensive  remodifications  to  allow  it  to  simulate  low 
altitude  orbits,  as  is  common  for  Earth  Observation  satellites. 

The  Soviet  density  model  gives  a  reasonable  first  estimate  of  the  confidence  in¬ 
terval  limits  expected  for  satellites  near  the  LANDSAT  6  altitude  of  700  km.  These 
limits  could  be  verified  and  improved  by  using  the  same  methodology  but  applying  a 
dynamic  atmospheric  model,  such  as  Jacchia-Roberts,  to  determine  the  atmospheric 
density  from  the  forecasted  and  actual  values  of  solax  flux  and  geomagnetic  activity. 
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Appendix  A 
Orbital  Elements 

A.l  Keplerian  Orbital  Elements 


The  Keplerian  orbital  element  set,  often  called  the  classical  orbital  elements,  are 
the  commonly  used  to  describe  an  orbit.  These  parameters,  stated  in  Table  A.l  allow 
the  orbit  to  be  ‘seen’  graphically.  Sometimes  the  true  anomaly,  /,  is  used  instead  of 
the  mean  anomaly,  M. 

Table  A.l.  Keplerian  Orbital  Elements 


a  =  semi- major  axis 

e  =  eccentricity 

i  =  inclination  of  the  orbit  with  respect  to  the 
equatorial  plane 

uj  =  argument  of  perigee  measured  from  the  line  of 
nodes  in  the  orbital  plane 

fl  =  longitude  of  the  ascending  node  measured  from  the 
vernal  equinox  in  the  equatorial  plane 
M  —  mean  anomaly  measured  from  perigee  to  mean 
satellite  position  in  orbit,  as  if  the  satellite 
had  constant  velocity  throughout  the  orbit  period 


The  semi-major  axis,  geometrically,  is  half  of  the  major  axis  of  the  elliptical  orbit. 
The  semi-major  axis  determines  the  size  of  the  orbit,  while  the  eccentricity  determines 
the  shape  of  the  orbit.  The  range  of  values  for  the  semi-major  axis  and  eccentricity 
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Table  A. 2.  Ranges  of  Semi-major  Axis  and  Eccentricity 


Conic  Type 

Semi-major  Axis 

Eccentricity 

Circle 

a  >  0 

e  =  0 

Ellipse 

a  >  0 

0  <  e  <  1 

Parabola 

a  =  oo 

e  =  1 

Hyperbola 

a  <  0 

e  >  1 

Table  A. 3.  Orbit  Types  from  Inclination 


Orbit  Type 

Inclination 

Polar 

i  =  90° 

Equatorial 

*  =  0  or  180° 

Direct 

i  <  90° 

Retrograde 

i  >  90° 

work  together  to  form  the  conic  shape  that  the  orbit  will  take.  These  ranges  for  the 
various  types  of  conic  orbits  are  shown  in  Table  A. 2. 

Perigee  is  the  point  closest  to  the  occupied  focus  of  the  elliptical  orbit,  ie.  the 
Earth.  The  line  of  nodes  is  defined  by  the  intersection  of  the  Earth’s  equatorial  plane 
and  the  satellite’s  orbital  plane.  It  is  at  this  point  that  the  angle  between  these  two 
planes  is  measured  to  determine  the  orbital  inclination.  The  inclination  is  always 
in  the  range:  0  <  i  <  180°.  To  avoid  any  ambiguity  of  which  angle  is  correct,  the 
inclination  is  measured  as  the  angle  between  the  perpendiculars  of  the  planes,  ie.  the 
z-axis  of  the  Earth  and  the  angular  momentum  vector  for  the  orbit  determined  by 
the  right  hand  rule.  Common  orbit  names  determined  by  the  type  of  inclination  are 
stated  in  Table  A. 3. 

The  vernal  equinox,  in  the  equatorial  plane,  is  always  defined  as  the  direction  to 
the  star  of  Aries.  The  Earth  Centered  Inertial  (ECI)  coordinate  system  uses  this  di¬ 
rection  to  orient  its  x-axis.  The  z-axis  is  defined  to  be  perpendicular  to  the  equatorial 
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plane  and  the  y-axis  is  defined  to  be  orthogonal  to  the  other  two  axes.  This  common 
coordinate  system  is  considered  sufficiently  inertial  for  satellite  systems  orbiting  the 
Earth.  An  illustration  of  the  ECI  coordinate  system  is  depicted  in  Figure  A.l. 


Z-axis 


Figure  A.l.  Earth  Centered  Inertial  Coordinate  System 

A. 2  Equinoctial  Orbital  Elements 

The  equinoctial  orbital  elements,  in  terms  of  the  classical  Keplerian  orbital 
elements  are  stated  in  Table  A. 4.  The  parameters  h  and  k  are  the  components  of 
the  eccentricity  vector  measured  from  the  vernal  equinox  in  the  orbit  plane.  The 
variable  A  is  the  ‘fast’  parameter  of  the  equinoctial  element  set  and  is  the  mean 
longitude  of  the  satellite.  The  equinoctial  orbital  elements  are  free  of  singularities 
in  the  Variation  of  Parameters  equations  of  motion.  The  common  Keplerian  element 
set  runs  into  division  by  zero  for  values  of  eccentricity  and  inclination  approaching 
zero  when  using  these  equations.  The  only  complication  for  the  equinoctial  element 
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Table  A. 4.  Equinoctial  Orbital  Elements 


set  exists  in  setting  a  retrograde  factor,  /.  This  retrograde  factor  can  equal  positive 
one  for  direct  equinoctial  elements  where  values  of  inclination  are  0°  <  i  <  180°. 
The  retrograde  factor  can  also  have  a  value  of  negative  one  for  retrograde  equinoctial 
elements  where  values  of  inclination  are  0°  <  i  <  180°.  The  direct  and  retrograde 
prefix  in  the  equinoctial  element  set  is  motivated  by  the  exact  limits  of  the  inclination. 
This  definition  is  different  than  the  Keplerian  element  set  where  the  general  range  of 
values  of  inclination  determines  the  prefix. 
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Appendix  B 


Porting  the  DSST  Standalone  to 
Non-IBM  Mainframe 
Environments 

B.l  VAX 

The  Draper  Semianalytic  Satellite  Theory  (DSST)  Standalone  is  stored  in  the  resi¬ 
dent  IBM  3090  mainframe  at  Draper  Laboratories  in  data  sets  LWE1 122. SAST. FORT 
and  LWEl  122. EVAL. FORT.  This  was  successfully  ported  to  a  VAX  8650  by  David 
Carter  [12j.  He  made  several  modulations  to  the  program  to  make  the  program 
compatible  with  VAX  FORTRAN. 

The  VAX  VMS  operating  system  is  not  compatible  with  file  names  with  the  char¬ 
acter#  in  tiiem.  For  greater  portability,  it  is  wise  not  tc  use  any  special  characters 
in  source  file  names. 

The  VAX  default  compiler  switch  for  double  precision  is  D -FLOATING.  This  is 
not  compatible  with  the  exponents  stated  in  the  SLP  data  sets,  to  avoid  this  problem, 
the  Standalone  must  be  compiled  using  the  /G-FLOATING  compiler  switch. 

In  a  CALL  statement,  hard  coded  constants  in  a  calling  sequence  can  not  be 
changed  inside  of  the  subroutine. 

Instead  of  using  DEFINE  FILE  to  define  a  unit  number  for  an  input  or  output 
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file,  the  OPEN  statement  is  the  correct  method  to  define  a  file  in  VAX  FORTRAN. 
Specific  to  the  DSST  Standalone,  units  14,  38,  and  78  are  the  only  direct  access  files 
needed  to  define  the  minimum  standalone  package. 

Quadruple  precision  float  variables  are  not  supported  in  standard  FORTRAN  77 
as  they  are  in  the  IBM  version  of  FORTR  AN.  In  the  Standalone,  these  variables  arc 
stated  but  not  used. 

Unique  to  the  Standalone,  two  subroutines,  RDNUMR  and  QREAD,  were  deleted 
since  no  other  subroutine  called  them  and  they  used  quadruple  precision  numbers. 

In  the  case  of  reading  a  direct  access  file,  specifically  in  subroutines  EVAL  and 
TCWF,  the  IBM  apostrophe  construction  must  be  replaced  by  ,REC=  . 

All  COMMON  and  EQUIVALENCE  declarations  must  precede  DATA  state¬ 
ments. 

All  BLOCK  DATAs  must  be  named  in  their  definition  in  every  subroutine. 

Specific  to  the  Standalone,  an  IMPLICIT  DOUBLE  PRECISION  (A-H,0- 
Z)  in  BLOCKDATA  MACHIN  was  needed  for  a  successful  run  of  the  Standalone. 

There  are  machine  specific  parameters  in  the  VAX  that  are  different  than  the  IBM 
mainframe,  ie.  minimum,  maximum,  null,  etc.  In  the  Standalone,  these  parameters 
are  in  BLOCKDATA  MACHIN  and  parameters  such  as  RELMAX,  RELMIN, 
and  RELNUL  need  to  be  changed.  The  values  specific  for  the  computer  that  is 
worked  on  can  be  found  in  its  manual  [12]. 

B.2  IBM  PC 

David  Carter  proved  that  it  is  a  very  simple  procedure  to  port  a  program  from  an 
IBM  mainframe  to  an  IBM  PC  using  Lahey  FORTRAN  [12],  To  support  this  port, 
the  following  listed  modifications  along  with  all  of  the  modifications  for  the  VAX  port 
were  needed. 
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For  proper  execution,  the  Lahey  FORTRAN  necessitates  the  use  of  FILE=...  in 
the  OPEN  statements. 

A  file  can  not  be  connected  in  Lahey  FORTRAN  on  the  IBM  PC  unless  it  has  first 
been  opened.  In  the  Standalone,  an  OPEN  and  CLOSE  statement  were  needed  in 
INREAD  to  connect  the  Newcomb  operator  file. 

If  the  Extended  Memory  version  of  Lahey  FORTRAN  is  used,  then  there  is  no 
need  to  downsize  the  COMMON  blocks.  This  version  of  the  Standalone  will  need  1.6 
Megabytes  to  run  successfully.  The  small  compiler  may  be  used  if  the  COMMON 
blocks  are  downsized  and  this  version  runs  with  only  600K  [12]. 

B.3  Sun  SPARCstation 

The  Draper  S^ir 'analytic  Satellite  Theory  (DSST)  Standalone  was  ported  from 
the  VAX  to  -  S  ■.  SPARCstation  1  (OS  version  4.03  Sun4C,  F77  Sun  release  4.0)  by 
David  O.iter  with  help  from  David  Kang.  This  port  took  approximately  2  hours. 
The  Standalone  default  test  case  output  from  the  Sun  was  compared  with  the  output 
from  the  VAX  and  the  data  agreed  to  the  twelfth  decimal  place.  The  modifications 
to  successfully  port  the  Standalone  are  listed  in  the  paragraphs  below. 

Any  RECL  specifications  in  direct  access  OPEN  statements  must  be  changed 
from  4-byte  convention  (required  by  VAX)  to  1-byte  convention  (standard). 

The  Sun  requires  that  FORTRAN  files  end  with  a  .f  rather  than  a  .for. 

The  data  files,  SLP1950,  SLPTOD,  NEWCOMB,  and  TIMECOEF,  must 
be  rewritten  so  that  they  are  in  the  correct  Sun  form.  This  is  accomplished  by 
running  the  rewrite  programs  writslp.f,  writtim.f,  and  writnukes.f.  On  the  Sun, 
the  SLP  program  can  be  compiled  and  linked  in  one  step  by  typing  f77  writslp.f, 
with  creates  an  executable  file  named  a. out  and  can  be  run  by  typing  its  name  and 
pressing  RETURN.  This  procedure  is  repeated  for  each  of  the  remaining  data  files 
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needed. 

Sun  FORTRAN  does  not  allow  a  block  data  subprogram  to  have  the  same  name  as 
its  associated  COMMON  block.  The  name  of  each  BLOCK  DATA  was  renamed 
to  end  in  the  character  $. 

All  of  the  subroutines  were  compiled  but  not  linked  by  using  the  wild  card  com¬ 
mand  f77  -c  *.f  after  temporarily  changing  the  name  of  the  main  program,  semi- 
anal. f  so  that  is  would  not  be  compiled  with  the  others. 

Large  subroutines  may  need  to  be  compiled  using  the  compiler  switch  -Nq999  if 
a  regular  compilation  fails  because  of  symbol  table  overflow.  This  compiler  option 
allows  more  symbols. 

The  main  program  was  compiled  and  linked  with  all  of  the  subroutines  using  the 
command  f77  *.o  semianal.f.  The  executable  code,  a.out  was  renamed  semianal 
with  the  command  mv  a.out  semianal.  The  output  from  a  run  was  commanded  to 
be  printed  to  a  file  named  aero. out  by  typing  semianal  aero. out. 

All  of  the  source  code,  object  files,  binary  data  files,  and  the  executable  code  was 
written  to  a  cartridge  tape  using  the  command  tar  -cv  *  [11]. 

B.4  Macintosh  PC 

There  are  many  FORTRAN  compilers  on  the  market  for  the  Macintosh  computer, 
but  the  only  one  found  that  could  manage  a  large  number  of  subroutines  was  Absoft’s 
MacFortran  II  that  is  used  as  a  tool  in  the  Macintosh  Programmers  Workshop  (MPW) 
from  Apple.  Others  that  failed  in  handling  a  large  program  ported  from  the  VAX 
were  Absoft’s  MacFortran  in  MPW  and  Mactran  Plus  from  DCM  Data  Products. 
Mactran  Plus  by  DCM  Data  Products  had  several  bugs  in  it  that  severely  limited 
the  number  of  subroutines  and  libraries  allowed.  In  addition,  their  Technical  Help 
Department  was  not  well  organized,  knowledgeable,  nor  understandable.  MacFortran, 
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by  Absoft,  which  also  runs  in  MPW  like  MacFortran  II,  had  bugs  that  interfered  with 
data  stored  in  common  blocks  [46], 

B.4.1  To  Install  MPW 

The  MPW  MacFortran  II  program  can  be  bought  from  Absoft  Corporation  for 
$595.00.  To  install  the  application  on  the  hard  drive  of  a  Macintosh  computer,  first 
copy  the  Installation  Folder  from  the  first  disk  to  the  folder  named  MPW,  on  the 
hard  drive. 

Open  the  Installation  Folder  on  the  hard  drive  and  launch  the  file  MPW  Installer 
(Launch  Me) 

Click  the  NO  button  at  the  Update  Alert  to  install  MPW. 

Insert  and  confirm  each  of  the  remaining  MPW  disks  (all  except  the  first  disk) 
and  click  the  YES  button  at  More  MPW  Disks  dialogs  that  will  appear.  Each 
file,  as  it  is  being  copied  to  the  hard  drive,  is  listed  and  each  disk  will  be  ejected 
automatically. 

Click  the  NO  button  at  the  More  MPW  Disks  dialog  after  installing  the  last 
of  the  MPW  disks. 

Click  the  OK  button  at  the  Update  Complete!  dialog. 

Throw  away  the  Installation  Folder  by  dragging  its  icon  to  the  trash  can.  It  is 
OK  to  throw  away  the  applications  in  the  Installation  Folder. 

Open  the  MPW  Folder  on  the  hard  drive  and  launch  the  file  MPW  Shell. 

Type  SetupMacFortran  on  a  line  by  itself  in  the  Worksheet  window  followed 
by  pressing  ENTER. 

At  the  Remove  Build  and  Project  menus?  dialog,  click  the  NO  button. 

Quit  MPW  by  selecting  the  QUIT  command  in  the  FILE  pull  down  menu  [47]. 


B.4.2  Changes  Needed  to  Port  a  Program 


To  port  a  program  from  a  VAX  to  the  Macintosh,  there  are  a  few  changes  that 
will  be  needed  to  get  it  running  in  MacFortran  II. 

All  of  the  common  blocks  must  be  the  same  size  as  their  associated  block  data 

set. 

Any  RECL  specifications  in  direct  access  OPEN  statements  must  be  changed 
from  Tbyte  convention  (required  by  VAX)  to  l-byte  convention  (standard)  as  in  the 
Sun  SPARCstation. 

The  reserved  unit  numbers  for  standard  input  and  output  are  listed  in  Table  B.l. 


Table  B.l.  MacFortran  II  Reserved  Unit  Numbers 


UNIT  NUMBER 

STANDARD 

5 

input 

6 

output 

9 

i/o 

* 

i/o 

It  is  easier  to  manage  a  large  program  that  will  be  repeatedly  modified  if  the 
compile  and  link  statements  are  kept  in  a  .make  file  that  can  be  rebuilt  any  time.  To 
do  this,  just  use  the  BUILD  pull  down  menu  and  select  MAKE  BUILD  COM¬ 
MANDS...  then  double  click  on  the  files  that  you  want  linked  together.  This  will 
automatically  create  the  .make  file  and  the  compiler  and  linker  options  will  have  to 
be  manually  added.  This  is  convenient  if  a  subroutine  is  modified,  since  the  program 
can  be  easily  rebuilt  by  selecting  BUILD  under  the  BUILD  pull  down  menu  and 
typing  in  the  name  of  the  .make  file  and  the  .make  file  will  automatically  recompile 
any  subroutines  that  have  been  modified  since  the  last  build  and  then  will  relink  them 
together  and  create  executable  code  or  a  program  application.  The  only  change  that 


134 


will  have  to  be  made  to  the  subroutines  to  accomplish  this  is  all  FORTRAN  files  must 
end  in  .f,  not  .for  as  in  the  Sun  SPARCstation. 

The  block  datas  in  separate  files  must  be  listed  first  in  the  link  line  listing  whether 
a  manual  link  or  a  link  in  a  .make  file  is  being  applied. 

The  compiler  options  used  to  successfully  port  a  program  from  the  VAX  are  listed 
in  Table  B.2.  Checking  the  array  boundaries  and  enabling  the  alignment  warnings  are 


Table  B.2.  Compiler  Options 


Option 

Description 

-C 

Check  array  boundaries 

-A 

Enable  alignment  warnings 

-N8 

Autosegmentation 

-f 

Fold  to  lower  case 

-s 

Static  storage 

not  necessary  but  are  beneficial  to  have  when  first  trying  to  port  the  program.  They 
ensure  that  the  disk  space  allocated  for  the  program  stays  intact  and  does  not  overflow 
into  other  applications  if  an  array  exceeds  its  boundaries.  The  autosegmentation  is 
needed  for  a  program  with  many  subroutines  of  varied  lengths.  If  this  option  is  not 
used,  the  linker  may  give  an  error  of  a  jump  step  being  too  large.  The  fold  to  lower 
case  and  static  storage  options  are  needed  by  programs  ported  from  the  VAX  and 
IBM  mainframe.  The  MacFortran  II  compiler  is  case  sensitive  to  the  type  in  variable 
names  unless  fold  to  lower  case  is  used  and  the  storage  of  variables  in  subroutines  is 
different  from  the  VAX  and  IBM  mainframe  if  the  static  storage  option  is  not  used. 

One  other  compiler  option  that  can  be  used  is  -N9.  This  option  allows  for  more 
checks  for  the  propellor  period  command  to  stop  executing  a  program.  This  option 
is  not  compatible  if  the  program  contains  computed  GOTO’s. 

The  linker  options  that  were  successfully  used  are  listed  in  Table  B.3. 
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Table  B.3.  Linker  Options 


Option 

Description 

-o  name 

Name  of  executable  file 

-d 

Suppress  duplicate  symbol  definition 

-f 

Largest  common  used 

-mf 

Use  Multifinder  memory 

-ss  1000000 

Segment  size  maximum 

-srt 

Sort  near /far  files 

After  the  program  works,  optimization  options  can  then  be  added.  These  opti¬ 
mizations  will  increase  compiler  time  but  will  decrease  runtime.  The  compiler  options 
-C  and  -A  can  be  deleted  to  decrease  the  compiler  time  after  the  program  works. 

Large  tools,  ie.  straight  executable  code,  made  by  MPW  share  memory  with 
the  MPW  application,  thus  they  may  not  even  be  able  to  load  -  so  it  is  better  to 
make  a  large  program  into  an  application.  Also,  the  default  for  applications  made  in 
MPW  is  that  they  are  automatically  linked  with  the  MacFortran  Runtime  Window 
Environment  (MRWE).  If  any  data  is  printed  to  the  screen,  a  window  is  automatically 
opened  with  save  and  print  menu  items  at  the  top  of  the  screen. 

A  PAUSE  rather  than  a  STOP  at  the  end  of  the  program  for  the  application  is 
helpful  if  any  data  is  written  to  the  screen.  By  adding  a  pause,  the  window  will  stay 
open  so  that  the  data  can  be  viewed  and  either  saved  in  a  file  or  printed  by  using  the 
menu  commands. 

If  an  application  is  not  running  correctly,  try  allocating  more  memory  to  the 
application  by  using  the  Get  Info  command. 

The  debugger  SADE  is  extra.  If  it  is  used,  SADE  can  only  run  under  Multifinder 
but  large  programs  can  usually  only  be  compiled  and  linked  under  Finder. 

Source  lines  can  be  shown  by  clicking  on  an  error  message  during  compilation. 
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f77  Option-X.f  compiles  all  .f  files  in  the  current  folder  and  creates  an  application. 

Option-X  matches  any  string  of  0  or  more  characters.  ?  matches  any  single 
character. 

Search  -i/  word  /  file...  searches  through  files  for  word  and  prints  out  the  file 
and  line  number. 

Option-d  is  the  line  continuation  character  for  MPW  commands  typed  into  the 
worksheet  [46]. 

B.4.3  Port  of  DSST  Standalone 

The  Draper  Semianalytic  Standalone  Orbit  Propagator  (DSST  Standalone)  was 
ported  from  the  IBM  3090  to  the  VAX  as  discussed  above.  Downloaded  from  the 
VAX,  the  Standalone  was  ported  to  a  Macintosh  Ilci.  The  DSST  Standalone  was 
successfully  run  using  the  default  test  case.  The  Macintosh  Ilci  output  was  manually 
checked  with  output  from  the  VAX  and  the  values  were  found  to  correlate  up  to 
the  ninth  decimal  place.  The  porting  was  finalized  on  December  11,  1990  after  six 
month’s  work  of  trial  and  error.  Now  that  a  FORTRAN  compiler  for  the  Macintosh 
has  been  found  to  work,  a  port  of  another  program  would  only  take  a  few  hours. 

The  DSST  Standalone  and  the  needed  input  files  can  be  stored  on  six  high  den¬ 
sity  disks.  The  Standalone  application  has  been  modified  to  include  namelist  inputs 
and  this  section  explains  how  to  install  the  application  in  a  Macintosh  and  use  it 
successfully. 

B.4.4  Needs  for  Exporting 

After  installing  MPW  version  3.1  and  MACFORTRAN  II  on  the  computer  (four 
disks),  create  a  new  folder  inside  of  the  MPW  folder  and  name  it  STANDALONE 
f.  Inside  of  this  folder  put: 
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The  five  name  list  files: 


NLSTELM  (7K) 
NLSTFRC  (7K) 
NLSTAVR  (7K) 
NLSTSP  (7K) 
NLSTDRV  (7K) 


The  four  input  data  files: 

NEWCOMB  (315K) 

SLP1950  (977K) 

SLPTOD  (977K) 

TIMECOEF  (21K) 

The  application  : 

STANDALONE  (343K) 

If  the  actual  FORTRAN  files  are  to  be  included  for  modification  by  the  user  -  be 
sure  to  include  the  build  file,  STANDALONE. make  (28K).  All  of  the  files  needed 
for  exporting  and  the  actual  FORTRAN  files  can  be  contained  on  six  disks. 

B.4.5  How  to  use  .make  File  to  Build  Application 

This  program  must  be  built  under  Finder  (as  opposed  to  Multifinder)  or  MPW 
will  not  be  able  to  run  the  application. 

Double  click  on  the  MPW  SHELL  icon  and  MPW  will  open  up  a  window  with 
menus  at  the  top. 

Pull  down  the  DIRECTORY  menu  and  select  SET  DIRECTORY. 

Click  once  on  the  folder  in  which  the  application  is  stored  (In  my  setup  this 
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is  STANDALONE  f)  and  then  click  once  on  the  DIRECTORY  button  (it  is 
highlighted). 

Pull  down  the  BUILD  menu  and  select  BUILD.  This  can  also  be  done  by  press¬ 
ing  Propellor-B. 

MPW  will  then  display  a  prompt  window  to  name  the  .make  file  to  be  built,  type 
STANDALONE  and  then  press  RETURN.  This  will  compile  any  subroutines  that 
have  been  changed  or  do  not  have  existing  object  files.  It  will  then  link  the  files  and 
create  the  application  STANDALONE. 

After  the  build  is  complete,  it  will  return  to  the  MPW  window  and  just  press 
ENTER  to  start  the  application  within  MPW. 

B.4.6  How  to  Run  Inside  MPW 

This  program  must  be  run  under  Finder  (as  opposed  to  Multifinder)  or  MPW 
will  not  be  able  to  run  the  application. 

Double  click  on  the  MPW  SHELL  icon  and  MPW  will  open  up  a  window  with 
menus  at  the  top. 

Pull  down  the  DIRECTORY  menu  and  select  SET  DIRECTORY. 

Click  once  on  the  folder  in  which  the  application  is  stored  (In  my  setup  this 
is  STANDALONE  f)  and  then  click  once  on  the  DIRECTORY  button  (it  is 
highlighted). 

Once  the  shell  returns  you  to  the  window  type  in  the  name  STANDALONE 
and  then  press  ENTER. 

The  information  will  be  printed  out  to  the  screen  and  you  may  send  this  to  the 
printer  at  the  end  of  the  program  by  pulling  down  the  FILE  menu  and  selecting 

PRINT  WINDOW. 

To  exit  after  the  program  is  finished  just  press  RETURN. 
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B.4.7  How  to  Run  as  a  Standalone  Application 

Once  again,  this  program  must  be  run  under  Finder  (as  opposed  to  Multifinder) 
or  this  program  can  not  be  opened. 

Double  click  on  the  STANDALONE  icon  and  MPW  will  automatically  start  up 
the  program. 

The  information  will  be  printed  out  to  the  screen  and  you  may  send  this  to  the 
printer  at  the  end  of  the  program  by  pulling  down  the  FILE  menu  and  selecting 

PRINT  WINDOW. 

To  exit  after  the  program  is  finished  just  press  RETURN. 

B.4.8  Input  Through  Namelists 

Any  data  that  needs  to  be  changed  from  the  DSST  standard  test  case  setting 
may  be  set  in  the  namelist  files.  A  namelist  can  be  opened  from  the  MPW  shell 
window  by  typing  open  nameofnamelist  and  then  pressing  ENTER.  Remember 
that  the  directory  must  first  be  set  to  the  folder  that  the  namelists  are  stored  in, 
ie.  STANDALONE  f.  A  namelist  file  could  also  be  opened  with  by  pulling  down 
FILE  menu  and  choosing  OPEN....  MPW  will  then  open  up  a  dialog  window  in 
which  the  file  to  be  opened  may  be  selected  and  opened  by  clicking  on  the  OPEND 
button  which  is  highlighted.  For  reference  to  variable  descriptions  of  the  variables 
that  may  be  set  in  the  namelist  files  see  Table  B.4.  One  common  variable  that  will  be 
set  is  the  epoch  orbital  elements.  These  are  set  in  namelist  NLSTELM.  The  default 
form  is  for  these  orbital  elements  to  be  Keplerian.  They  are  stored  in  an  array  in  the 
following  order,  semi-major  axis,  eccentricity,  inclination,  longitude  of  the  ascending 
node,  argument  of  perigee,  and  mean  anomaly. 

The  namelist  input  format  is  of  the  following  form: 

$name  of  namelist 
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Table  B.4.  Variable  Reference 


NAMELIST  FILE 

VARIABLE 

DESCRIPTION 

NLSTELM 

NLSTAVR 

NLSTFRC 

NLSTDRV 

NLSTSP 

setelm.f 

setavr.f 

setfrc.f 

setdrv.f 

setsp.f 

member  =  value, 
member  =  value, 
member  =  value, 

SEND 

The  end  of  the  data  must  be  followed  by  a  carriage  return  so  that  the  end  of  the 
file  may  be  detected  when  executing. 

B.4.9  After  Application  is  Finished 

After  the  application  is  finished,  there  is  a  pause  at  the  end  so  RETURN  must 
be  pressed  to  go  back  to  the  MPW  Shell  window.  Before  RETURN  is  pressed,  this 
is  the  only  chance  to  view  the  output  from  the  execution  and  either  save  and/or  print 
it. 

To  quit  the  MPW  Shell,  pull  down  the  FILE  menu  and  select  QUIT  [47]. 
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